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Abstract. We discuss the dynamics and thermodynamics of the Hamiltonian Mean Field model (HMF) 
which is a prototypical system with long-range interactions. The HMF model can be seen as the one Fourier 
component of a one- dimensional self-gravitating system. Interestingly, it exhibits many features of real self- 
gravitating systems (violent relaxation, persistence of metaequilibrium states, slow collisional dynamics, 
phase transitions,...) while avoiding complicated problems posed by the singularity of the gravitational 
potential at short distances and by the absence of a large-scale confinement. We stress the deep analogy 
between the HMF model and self-gravitating systems by developing a complete parallel between these two 
systems. This allows us to apply many technics introduced in plasma physics and astrophysics to a new 
problem and to see how the results depend on the dimension of space and on the form of the potential of 
interaction. This comparative study brings new light in the statistical mechanics of self-gravitating systems. 
We also mention simple astrophysical applications of the HMF model in relation with the formation of 
bars in spiral galaxies. 

PACS. 05.20.-y Classical statistical mechanics - 05. 45. -a Nonlinear dynamics and nonlinear dynamical 
systems 



1 Introduction 

The statistical mechanics of systems with long-range in- 
' teractions is currently a topic of active research in physics 
\ because it differs in many respects from that of more fa- 
. miliar systems with short-range forces that are extensive 
' (Dauxois et al. 2002a). Among long-range interactions, 
, gravity is probably the most important and most funda- 

■ mental example (Padmanabhan 1990, Chavanis 2002a). 
\ However, the statistical mechanics of self-gravitating sys- 

tems initiated by Antonov (1962) and Lynden-Bell (1968) 

■ is complicated due to the divergence of the gravitational 
, force at short distances and to the absence of shielding 

(or confinement) at large distances. These difhculties are 
specific to the gravitational force and not to the long- 
range nature of the interaction. Therefore, it may be of 
conceptual interest to consider simpler systems with long- 
range interactions to distinguish what is specific to the 
gravitational force and what is common to systems with 
long-range interactions. 

A toy model of systems with long-range interactions 
is the so-called HMF (Hamiltonian Mean Field) model. It 
consists of N particles moving on a circle and interacting 
via a cosine binary potential. This can be seen as a one- 
dimensional plasma where the potential of interaction is 



truncated to one mode. This model is of great concep- 
tual interest because it exhibits many features present in 
more realistic systems with long-range interactions such as 
gravitational systems. In addition, it is sufficiently simple 
to allow for accurate numerical simulations and analytical 
results. 

To our knowledge, what is now called the HMF model 
was first introduced by Konishi & Kaneko (1992). They 
found that a cluster is formed in some cases and that the 
system remains uniform in other cases. Inagaki & Konishi 
(1993) realized that the Konishi-Kaneko system is nothing 
but the one Fourier component of a one-dimensional self- 
gravitating system and explained the formation of clusters 
as an instability similar to the Jeans instability in self- 
gravitating systems described by the Vlasov equation. In- 
agaki (1993) studied the thermodynamical stability of the 
Konishi-Kaneko system and identified the existence of a 
second order phase transition at a critical temperature T^. 
Above Tc the only statistical equilibrium state is uniform, 
whereas below Tc this uniform state looses its thermody- 
namical stability and clustered states appear. In order to 
justify his results dynamically, Inagaki (1996) developed a 
"collisional" kinetic theory of the Konishi-Kaneko system 
based on results of plasma physics and proposed to model 
the dynamics of the system by the Lenard-Balescu equa- 
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tion for a one dimensional plasma truncated to one mode. 
However, as we shall see, his conclusions demand further 
discussion. 

The same model was considered at about the same 
time by Pichon and Lynden-Bell (Pichon 1994) who gave 
an astrophysical application of this model in relation with 
the formation of bars in galactic disks. In their approach, 
the stars follow rigid elliptical orbits with eccentricity e. If 
represents the inclination of ellipse i and f2i its angular 
velocity, the torque exerted by an orbit to the other can 
be written a~^dflx/dt = d^i2/d(f)i where is the adi- 
abatic moment of inertia of the inner Lindblad orbit and 
?/'i2 = GA'^ cos 2(01 — (/)2) is the effective alignement poten- 
tial. At high temperatures, the orbits are almost uniformly 
distributed in space and the system is in a disk phase (see 
Fig. 01). However, below some critical temperature Tj., the 
ellipses tend to align to each other and form a bar (see 
Fig. 13). Those bars are reported observationally in real 
galaxies. Pichon and Lynden-Bell studied the linear sta- 
bility of these bars with respect to the Vlasov equation 
and proposed that the clustered phase could result from 
a process of violent relaxation, a concept introduced by 
Lynden-Bell (1967) to explain the regularity of collision- 
less stellar systems such as elliptical galaxies. 

The Konishi-Kaneko system, now called the HMF model, 
also appeared in statistical mechanics (Antoni & Ruffo 
1995). In that context, the motivation was to devise a sim- 
ple model with long-range interactions keeping the rich- 
ness of more realistic systems while being amenable to a 
full analytical and numerical treatement. Excitingly, this 
simple model displays a lot of interesting features (vio- 
lent relaxation, persistence of metaequilibrium states, slow 
collisional relaxation, phase transitions,...) also present in 
other systems with long-range interactions such as stellar 
systems and 2D vortices (Chavanis 2002a). The proper- 
ties of the HMF model have been studied in great detail 
in a lot of recent papers (see Dauxois et al. 2002b for a re- 
view). Despite its oversimplification, the HMF model can 
be seen as a pedagogical model to take a step into the 
physics of systems with long-range interactions. It is said 
sometimes to represent the "harmonic oscillator" of sys- 
tems with long-range interactions. This probably explains 
its popularity. 

In the present paper, we shall emphasize the connec- 
tion between the HMF model and the results established 
in astrophysics and plasma physics. In particular, we will 
adapt the methods developed for 3D self-gravitating sys- 
tems to the case of a one-dimensional system of particles 
with cosine interactions. The motivation of this extension 
is two-fold. The first is to show that the results obtained 
in astrophysics and plasma physics can have applications 
in other domains of physics, including the HMF model 
(this has not been sufficiently appreciated by workers in 
that field since the early work of Inagaki) . The second is 
to stress the analogies and the differences which appear in 
long-range systems as we change the dimension of space 
and the potential of interaction. Among the analogies be- 
tween 3D self-gravitating systems and the HMF model, 
we note: the concept of violent relaxation and the slow 



collisional dynamics. Among the differences, we note: the 
equivalence of statistical ensembles for the HMF model 
(contrary to 3D gravitational systems), the existence of 
second order phase transitions (instead of first order or ze- 
roth order phase transitions for 3D gravitational systems) 
and the vanishing of the collision operator at the order 
1/7V in the BBGKY hierarchy contrary to the Coulom- 
bian or Newtonian case. 

The paper is organized as follows. In Sec. 12 we con- 
sider the statistical equilibrium states of the HMF model 
in both microcanonical and canonical ensembles. We syn- 
thesize previously known results and we derive explicit cri- 
teria of thermodynamical stability for the uniform phase 
as well as for the clustered phase. We also describe correc- 
tions to the mean-field approximation close to the critical 
point. In Sec. we consider a one-dimensional gaseous 
system with cosine interactions (the analogue of a "gaseous 
star" ) described by the Euler equations with a barotropic 
equation of state. We discuss in particular the equivalent 
of the Jeans instability. In Sec. 01 we consider the coUi- 
sionless evolution of the HMF model (the analogue of a 
"stellar system") described by the Vlasov equation and 
discuss the concept of violent relaxation and metaequilib- 
rium states. We interprete these quasi-equilibrium states 
as particular stationary solutions of the Vlasov equation 
on the coarse-grained scale resulting from phase mixing 
and incomplete violent relaxation. We regard Tsallis func- 
tional Sq[f] = — ^ / {f^ — f)dddv and Boltzmann func- 
tional Ssif] = — / f\nfd9dv as particular _ff-functions in 
the sense of Tremaine et al. (1986) associated with poly- 
tropic and isothermal distributions. We study the dynami- 
cal stability of stationary solutions of the Vlasov equation 
and compare with the dynamical stability of stationary 
solutions of the barotropic Euler equations. This is the 
same type of comparison as between "gaseous systems" 
and "stellar systems" in astrophysics. In that respect, we 
discuss the equivalent of the Antonov first law (Binney & 
Tremaine 1987) for the HMF model. We derive a crite- 
rion of nonlinear dynamical stability for steady states of 
the Vlasov equation of the form / = /(e) with /'(e) < 
where e is the individual energy, and show that it can 
be written as a condition on the velocity of sound in the 
corresponding barotropic gas. This criterion is equivalent 
to the criterion obtained by Yamaguchi et al. (2004) but 
it is expressed differently. We also analyze the linear dy- 
namical stability of steady states of the Vlasov equation 
and study the dispersion relation for isothermal and poly- 
tropic distributions. In Sec. El we discuss the collisional 
evolution of the HMF model and explain why the ki- 
netic theory is more complicated than for 3D Newtonian 
interactions. In particular, the Landau and the Lenard- 
Balescu collision terms vanish for ID systems so that the 
evolution of the system as a whole is due to terms of 
order smaller than 1 /N in the expansion of the correla- 
tion functions for N — )■ +oo. This implies that the relax- 
ation time is larger than Ntn (where tu is the dynam- 
ical time). By contrast, we can develop a kinetic theory 
at order 1/A^ to analyze the relaxation of a "test parti- 
cle" in a bath of "field particles" with static distribution 
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/o(u). The evolution of the distribution function P{v, t) of 
the velocity of the test particle satisfies a Fokkcr-Planck 
equation. We give explicit expressions for the diffusion co- 
efficient and the auto-correlation function and we com- 
pare their expressions depending on whether collective ef- 
fects are taken into account or not. We also show that the 
auto-correlation function decreases exponentially rapidly 
in time with a rate coinciding with the damping rate 7 of a 
stable perturbed solution of the Vlasov equation. Finally, 
in Sec. El we discuss the case of self- attracting Brownian 
particles described by non-local Fokker-Planck equations. 
This stochastic model is the canonical counterpart of the 
Hamiltonian TV-body problem. It could be called the BMF 
(Brownian Mean Field) model. We study the dynamical 
stability of steady states of the non-local Smoluchowski 
equation and solve this equation numerically to show the 
formation of a clustered state from an unstable homoge- 
neous state due to long-range interactions. We also provide 
analytical solutions of the dynamics close to the critical 
point Tc and for T = 0. All these models have an equiva- 
lent in the astrophysical literature and this paper stresses 
the analogies and differences between self-gravitating sys- 
tems and the HMF model. This comparative study brings 
new light in the statistical mechanics of self-gravitating 
systems by showing what is specific to gravity and what 
is common to systems with long-range interactions. 



2 Statistical equilibrium 

2.1 The mean-field approximation 

We consider a system of N particles moving on a circle 
and interacting via a cosine binary potential. This is the 
so-called HMF model. As explained in the Introduction, 
this model can also describe a system of stars moving on 
elliptical orbits, each orbit exerting a torque on the others. 
Fundamentally, the dynamics of this system is governed 
by the Hamilton equations 



d0, _ dH 
' dt dvi ' 



' dt 



dH 



N 1 k 

1=1 

where 9i is the angle that makes particle/ellipse i with an 
axis of reference and k is the coupling constant (similar 
to the gravitational constant G). In the rest of the paper, 
we shall refer to this system as a "stellar system" ; this is 
to emphasize the analogies with real 3D stellar systems 
whose dynamics is also governed by Hamiltonian equa- 
tions with long-range interactions. We have also general- 
ized the usual HMF model to a population of particles 
with different masses m^. However, in most of the paper, 
we shall assume that all the particles have the same mass 
m — 1. The multi-species HMF model will be discussed 
specifically in Sec. [3 



The evolution of the A^-body distribution function is 
governed by the Liouville equation 



dP, 



N 



dt 



+ 



N 

E 

=1 



dP, 



N 



dP, 



N 



dvi 



= 



(2) 



where Fi = —■^^f^iSm{9i — 6j) is the force experi- 
enced by particle i. Any distribution of the form Pn — 
x{H)5{E — H) is a stationary solution of the Liouville 
equation. For A^ 1 (fixed) and t +00, this system 
is expected to reach a statistical equilibrium state due to 
the developement of correlations between particles (this 
will be refered to as a "collisional" relaxation). As is cus- 
tomary in statistical mechanics, we shall assume that the 
equilibrium A^-body distribution function is described by 
the microcanonical distribution 



Pat (6*1, wi, 



,vn) = 



1 



9{E) 



6{E- H), 



(3) 



expressing that all accessible microstates (with the right 
values of energy and mass) are equiprobable. Whether this 
is indeed the case has not been proved rigorously as it 
relies on a hypothesis of ergodicity, so this statement is 
essentially a postulate. 

For systems with long-range interactions (self-gravitating 
systems, 2D vortices, HMF model,...), it can be shown 
that the mean-field approximation is exact in an appro- 
priate thermodynamic limit. This can be shown for ex- 
ample by considering an equilibrium BBGKY-like hier- 
archy (Chavanis 2004b,2005). For the HMF model, the 
thermodynamic limit is A^ — > -|-cxd in such a way that 
the properly normalized energy e = SttE /kM'^ and tem- 
perature rj = pkM/An are fixed, where M = Nm is the 
total mass. These control parameters are similar to those, 
e = ER/GM^ and tj = l3GMm/R, describing 3D gravi- 
tational systems (see, e.g., Chavanis 2003a). In that limit 
A^ +00, the two-body distribution function can be ex- 
pressed as a product of two one-body distribution func- 
tions 

P2iOi,Vu02,V2) = Pliei,Vi)Piie2,V2) + Oil/N). (4) 

The average density of particles in phase space is given by 
/(^,«) = - e,)S{v - «,)) = NPi{0,v). The total 

mass can then be expressed as 

M = J fdOdv. (5) 

On the other hand, the average energy E = {H) is 

E ^ N Pi{e,v)—dedv 

--^N{N~l) I cos{e-e')P2i9,v,e',v')d9dvd9'dv'.{6) 
in J 

In the mean-field limit, it reduces to 



f fv^dOdv + i / f-^dOdv, 



(7) 
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where 

m = ^^l'^ cos{0-e')p{0')de', (8) 

is the potential and p = J fdv is the spatial density. Note 
that the average force experienced by a particle located in 
6 is {F) = —<P'{0). If p is symmetric with respect to the 
X-axis, so that p{^0) = p{d), the foregoing relation can be 
rewritten 

^{e)^Bcos0, (9) 

where 

B = / p{0')cos0'd0'. (10) 

This parameter B is the equivalent of the magnetization 
(usually denoted M) in the case of spin systems. Inserting 
the relation ^ in Eq. ||7J), we find that the energy can be 
rewritten 

E=IJ fv'd0dv~^, (11) 

so that the potential energy is directly expressed in terms 
of B. 



2.2 The Boltzmann entropy 

We wish to determine the macroscopic distribution of par- 
ticles at statistical equilibrium, assuming that all accessi- 
ble microstates (with given E and M) are equiprobable. 
To that purpose, we divide the /i-space {9, v} into a very 
large number of microcells with size h. We do not put 
any exclusion, so that a microcell can be occupied by an 
arbitrary number of particles. We shall now group these 
microcells into macrocells each of which contains many mi- 
crocells but remains nevertheless small compared to the 
phase-space extension of the whole system. We call ly the 
number of microcells in a macrocell. Consider the configu- 
ration {ui} where there are ni particles in the I''' macro- 
cell, 712 in the 2"*^ macrocell etc. Using the standard combi- 
natorial procedure introduced by Boltzmann, the number 
of microstates corresponding to the macrostate {n^}, i.e. 
its probability, is given by 

W{{n,}) = ml[—. (12) 

i 

This is the Maxwell-Boltzmann statistics. As is customary, 
we define the entropy of the state {rii} by 

S{{n,})=\nW{{n,}). (13) 

It is convenient here to return to a representation in terms 
of the distribution function giving the phase-space density 
in the i-th macrocell fi — f{0i,Vi) = rii/vh. Using the 
Stirling formula and passing to the continuum limit v ^ Q, 
we obtain the usual expression of the Boltzmann entropy 

SbU]^- f flnfd0dv, (14) 



up to some unimportant additive constant. Then, the sta- 
tistical equilibrium state, corresponding to the most prob- 
able distribution of particles, is obtained by maximizing 
the Boltzmann entropy (|14|l at fixed mass M and energy 
E, i.e. 

Max {SbU] I E[f] = E, M[f] = A/}. (15) 

This maximization problem defines the microcanonical equi- 
librium state, which is the correct description for an iso- 
lated Hamiltonian system. 

We shall also consider the canonical description which 
applies to a system in contact with a thermostat impos- 
ing its temperature T. We will give an example of canon- 
ical system in Sec. corresponding to Brownian particles 
in interaction described by stochastic (not Hamiltonian) 
equations (BMF model). In the canonical ensemble, the 
statistical equilibrium state is obtained by minimizing the 
free energy ^ E[f] - TSbU] at fixed mass M and 

temperature T, i.e. 

Min {FbU] \M[f]^M}. (16) 

The relation between the Boltzmann entropy S'b[/] and 
the density of states g{E) in the microcanonical ensemble 
and between the Boltzmann free energy Fb [/] and the par- 
tition function Z(f3) in the canonical ensemble is discussed 
in Chavanis (2004b, 2005). The variational problems (fT5|l 
and (|16|l correspond to a saddle point approximation in 
the functional integral representation of g{E) and Z{f3). 

The variational problem (|15|l has been first investi- 
gated, in the HMF context, by Inagaki (1993). We review 
and precise the main results of his study and present an 
alternative derivation of the condition of thermodynami- 
cal stability using methods similar to those introduced by 
Padmanabhan (1990) and Chavanis (2002b,2003a) for 3D 
self-gravitating systems. This will make the analogy be- 
tween the two systems (stellar systems and HMF model) 
closer. This will also allow us to study the thermodynam- 
ical stability of the clustered phase, while the analysis of 
Inagaki (1993) is restricted to the uniform phase. 

2.3 The first variations: the Maxwell-Boltzmann 
distribution 

We need first to determine the critical points of entropy at 
fixed mass and energy. We write the variational principle 
in the form 

SSb - PSE - aSM = 0, (17) 

where f3 — I /T (inverse temperature) and a (chemical po- 
tential) are Lagrange multipliers enforcing the constraints 
on E and M . The solution of H17|l is the mean-field Maxwell- 
Boltzmann distribution 

/ = A'e-''(^+*), (18) 

where <!> depends on / through Eq. © . The foregoing rela- 
tion is therefore an integro-differential equation. Integrat- 
ing over the velocity, we obtain the mean-field Boltzmann 
distribution 

p = ^e-/3*. (19) 
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The same distributions p8(l and 119|l are obtained in the 
canonical ensemble by cancelling the first variations of free 
energy at fixed mass, using 6Fb ~ a5M = 0. Therefore, 
the critical points of the variational problems (|15|l and 
coincide. In the microcanonical ensemble, we need to 
relate the Lagrange multiplier (3 to the energy E. This 
defines a series of equilibria (3 = P{E). In the canonical 
ensemble, the inverse temperature (3 is assumed given and 
the corresponding (average) energy is obtained by invers- 
ing the graph l3{E). Of course, only the stable part of the 
series of equilibria is of physical interest, and defines the 
caloric curve (see Sec. 12 .511 . We note that the equilibrium 
distributions 1)18(1 and 1(19(1 can also be obtained from an 
equilibrium BBGKY-like hierarchy in the thermodynamic 
hmit N +00 defined in Sec.lOfChavanis 2004b,2005). 

Using the relation 0, the distribution of particles at 
statistical equilibrium is given by 



(20) 



The axis of symmetry is determined by the initial condi- 
tions. If i? = 0, the density p is uniform. This defines the 
homogeneous phase. If _B ^ 0, we have inhomogeneous 
states with one cluster at 6* = (if B < 0) or at = tt (if 
B > 0). The constant A is related to the mass by 



M = 2ttAIq{(3B), 
where /„ are the modified Bessel functions 



Iniz) = - 
TT 



' cos{n9)d9. 



For z ^ 0, 



Iniz) = {-zr 



r(7i + i) 4r(n + 2) 



+ 



and for 



In{z) 



1 - 



8z 



(21) 



(22) 



(23) 



(24) 



Using Eqs. I(1U() and ((22(1 we find that the order pa- 
rameter B is determined as a function of the temperature 
P by the implicit equation 



B 



kM Ii{/3B) 
27r Io{/3B)- 



(25) 



Setting X — (3B, we can rewrite the foregoing relation in 
the form 

(26) 



kM 



Then x, and consequently B, is determined as a function 
of T by a simple graphical construction sketched in Fig. 
nj We see that B = is always solution although B 
is possible only if 



T < 



kM 
An 



(27) 
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T>T„ 
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Fig. 1. Graphical construction showing the appearance of a 
clustered phase below some critical temperature Tc. 
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Fig. 2. Order parameter B (magnetization) as a function of 
the inverse temperature. 



In terms of the energy pi(l this corresponds to 



E < 



kM^ 
"8^ 



Er 



(28) 



The function B{T) is shown in Fig. [3 and its asymptotic 
behaviours are given in Sec. 12.41 Figure [3 displays a sec- 
ond order phase transition. We have a situation similar 
to a gravitational collapse below a critical temperature Tc 
or below a critical energy Ec- For T > T^,, the system is 
homogeneous. For T < Tc, the system forms one cluster 
around 9 = (for i? < 0) or around 9 = tt (for S > 0). At 
T = 0, the equilibrium state is a Dirac peak p = M6{9—7r) 
(for B = B,nax)- Density profiles are plotted in Fig. |31for 
different values oi x = [3B{0). Using the stellar disk inter- 
pretation of Pichon & Lynden-Bell, we have represented 
some stellar orbits in Figs. 01 and [3 by randomly choosing 
the orbits' angles with the equilibrium distribution p{6). 
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Fig. 6. Specific heat C = dE/dT as a function of temperature. 
Fig. 3. Evolution of the density profile as temperature is de- it experiences a discontinuity at the critical temperature T^. 
creased (from bottom to top). 
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Fig. 4. Stellar orbits in the "disk phase" for T > Tc. 
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Fig. 5. Stellar orbits in the "bar phase" for T < Tc. 



The "disk phase" for T > Tc is represented in Fig. ^ and 
the "bar phase" for T < Tc is represented in Fig. [S] 



2.4 The thermodynamical parameters 

According to Eq. H2ti|) . the relation between the temper- 
ature and the order parameter can be written in dimen- 
sionless form as 



2 Ii(x) 



For a: ^ 0, 



and for x +oo, 



X f 1 



Returning to original variables, we deduce that 
B ' 



(29) 



(30) 



(31) 



= ±W2(1-^ 



B 



(0<%^«l), (32) 



Tr. 



Br, 



T 



(33) 



where Bmax — "^Tc — kM/2'K is the maximum value of the 
magnetization obtained for T = when all the particles 
are at 9 = n. With this notation, the parameter x can be 
written 

X = 2TjB/B,nax- (34) 

On the other hand, for the Maxwellian velocity distri- 
bution p8(l . the expression of the energy 1)11(1 becomes 



E='-MT-^-^. 
2 k 



(35) 
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In terms of dimensionless parameters, we get 

1/ x'^' 
e = E E, = - 1 - - 
77 V 2?7 

For the homogeneous phase _B = 0, we simply have 

1 

e = — . 

■n 



(36) 



(37) 



For the inhomogeneous phase, we can easily obtain asymp- 
totic expansions. For a; — > 0, 



and for x 



-00, 



-2(1-- 

X 



Returning to original variables, we deduce that 

E _IT 
Eq 2Tc 



{T > n), 



Elf T 



E _ T ^ 

Eq Tc 



(0<^^«1), 
(T « T,), 



(38) 
(39) 

(40) 

(41) 
(42) 




£=87iE/kM 

Fig. 7. Calorique curve (series of equilibria) for the HMF 
model. The system displays a second order phase transition 
at a critical point Ec,Tc. 



E = ^MT + i / p<Pde. 



(46) 



We now take the variations of entropy around an equilib- 
rium solution. To second order 

2 



2 T 



T 



Sp\npd6 — — 



where —Eq — —kM"^ /Aw is the minimum value of energy 
obtained for T = (eo = -2). For T ^ T' , the specific 

heat C = dE/dT is given by C = |M and for T > by Now, the conservation of energy implies 



(Spy 



c 



M 



Therefore, at the critical point, it experiences a 



discontinuity (see Fig.lHI): 

C{T~) - C(T+) = 2M. 



Q = 5E= -M6T 
2 



<P5pde + - 



de. 

(47) 

(48) 



(43) 



The caloric curve/series of equilibria [3{E) is shown in Fig. 
[7| It displays a second order phase transition at (e, r?) = 
(1,1). This is different from 3D gravitational systems which 
rather display first order and zeroth order phase transi- 
tions (see, e.g., Chavanis 2002c). 



2.5 The second variations: thermodynamical stability 

To analyze the thermodynamical stability of the solutions 
determined by the variational problems (|15|l and (|16|l . we 
use an approach similar to that followed by Padmanab- 
han (1990) and Chavanis (2002b,2003a) in the case of 3D 
self-gravitating systems. We first maximize 5*3 [/] at fixed 
M[/], E[f] and p[e). This gives the Maxwelhan 



Eliminating 5T , we find that 
1 

~ ^2r 

1 



5^St 



Spd<Pd0 

\ 2 



<PSpd0 



{5pf 



de. 



We define the quantity q by the relation 

dq 



Sp^ 



d0' 



(49) 



(50) 



Physically, q — Spdd represents the mass perturbation 
within the interval [0,6]. Then, the conservation of mass 
is equivalent to q{0) = g(27r) = 0. Inserting this relation in 
Eq. 149(1 and using integrations by parts, we can put the 
second order variations of entropy in the quadratic form 



1 



P{0)e 



(44) 



Then, we can re-express the entropy and the energy as a 
function of the density in such a way that 



Sb = ^Af InT- / p\np d9, 



(45) 



with 

k 



dede'q{e)K{e,e')q{6'), (51) 



AtiT 



■ sm 



1 d<!> 
-9')— + -6(9-0') — 



dO' 



1 



p{6') de' 



(52) 
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Fig. 8. Dependence of the largest eigenvalue A with the 
temperature. A negative value of A corresponds to stability 
((5^5 < 0) and a positive value of A corresponds to instability. 



We are thus led to consider the eigenvalue problem 



27r 



(53) 



This yields 



d (\dq 
d9\J)d9 



27rT 



2ir 



q{9')cos{9-e')d6' 



2V d$ 



where 



do^^'^''- 



+ 2Ag, (54) 



(55) 



The system is stable if all A < and unstable if at least 
one A > 0. So far, we have worked in the microcanonical 
ensemble. If we work in the canonical ensemble, we have 
to minimize the free energy Fb = E — TSb at fixed mass 
and temperature. We can easily check that fixing T in the 
preceding calculations amounts to taking ^ = 0. Thus, 
instead of Eq. H54I) , we obtain 



d /I dq 
dd\J)d9 



2ttT 



9') cos(6' - e')d9' = 2\q. (56) 



2.6 The condition of thermodynamical stability 

If we consider the stability of the uniform solution p ■ 
M/2tt and = 0, the foregoing equations simplify into 



k 



27r d'^q 

Mde^ 2ttT 



2tt 



q{6') cos{e - 9')d9' = 2\q. (57) 



The eigenvalue equation is the same in the two ensembles. 
Hence, the stability criteria coincide, implying that the 
statistical ensembles (microcanonical and canonical) are 



equivalent. This is at variance with the case of 3D stellar 
systems (Padmanabhan 1990, Chavanis 2002b, 2003a). 

We can study the solutions of Eq. H57() by decompos- 
ing q in Fourier series. For the mode n, we have g„ = 
A,iSin(n0). For n 7^ 1, we get A„ = < showing 

that these modes do not induce instability. For n = 1, 
we have Ai = — -§[■ The uniform solution will be un- 
stable if Ai > yielding condition (|27|l . Therefore, the 
uniform phase is stable for T > Tc while it is unstable for 
T < Tc- By using the theory of linear series of equilib- 
ria (Katz 1978, 1980, 2003), applied here to a bifurcation 
point, we directly conclude from the inspection of Fig. 
that the clustered phase will be stable for T < Tc when 
the homogeneous phase becomes unstable. 

More precisely, it is possible to solve Eqs. l(5^ and 
(|56|l analytically for the clustered phase in the limit B — > 
0, which is valid close to the critical point {Ec,Tc). The 
calculations are detailed in Appendix A. In the canonical 
ensemble {V — 0), it is found that the largest eigenvalue 
is 



^27r 



f3c 



1 



and in the microcanonical ensemble (V 7^ 0) that 



XM 



(58) 



(59) 



More generally, the exact value of A obtained by solving 
Eqs. H54I) and H56II numerically is plotted versus the inverse 
temperature in Fig.|S| Since A < 0, we check explicitly that 
the clustered phase is stable. 



2.7 Correction to the mean-field approximation close 
to the critical point 

We can obtain the expression of the two-points correla- 
tion function from an equilibrium BBGKY-like hierarchy 
by closing the second equation of the hierarchy with the 
Kirkwood approximation (Chavanis 2004b, 2005). This is 
valid to order l/N in the thermodynamic limit defined 
previously. For the HMF model, it is then possible to ob- 
tain an explicit expression of the correlation function in 
the homogeneous phase. Writing the two-body distribu- 
tion function as 



it is found that 



KOi - 62) 



N 1 - P/Pc 



(60) 



(61) 



We note that the correlation function diverges close to the 
critical point (3 ^ 13c where the clustered phase appears 
and the homogeneous phase becomes unstable. This im- 
plies that the mean-field approximation ceases to be valid 
close to the critical point. We expect a similar result for 
3D self-gravitating systems although the situation is more 
difficult to analyze as (real) self-gravitating systems are 
always inhomogeneous. 
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If we take into account the contribution of non-trivial 
pair correlations (|61(l in the potential energy, we find fur- 
thermore that Eq. H37|l is replaced by 



1 2t] 



N 1 



(62) 



Therefore, finite N effects modify the shape of the caloric 
curve in the vicinity of the critical point. The mean- field 
approximation is valid if A^(l — 77) ^ 1. When the mean 
field approximation is valid, its order one correction for 
the specific heat is 



N 



1 



1 



1 



ttN {T/Tc - 1)2 



(T > Tc). (63) 



Finally, it is found that the spatial correlations of the 
force are given by (Chavanis 2004b, 2005) 



(mm) = 



1 



47r 1 - P/(3, 



C0S( 



In particular, the variance of the force is 



pfc2 



47r 1 - P/I3c 



(64) 



(65) 



Note that without the correlations, we would have sim- 
ply obtained (F^) ~ ^ which corresponds to the high 
temperature limit (T — > -l-oo) of Eq. H65() . 

For the HMF model, the variance (j65|l of the force is 
finite while the variance of the Newtonian force for 3D 
self- gravitating systems is infinite (Chandrasekhar & von 
Neumann 1942). For the HMF model, the distribution of 
the force is normal (Gaussian) while the distribution of 
the gravitational force in Z) = 3 is a particular Levy law 
called the Holtzmark distribution. On the other hand, for 
2D point vortices, the variance of the velocity is a marginal 
Gaussian distribution intermediate between normal and 
Levy laws (Chavanis & Sire 2000). Therefore, these three 
systems with long-range interactions (self-gravitating sys- 
tems, 2D vortices and HMF model) have their own speci- 
ficities despite their overall analogies. 



3 Gaseous systems 

As indicated in the Introduction, the HMF model is simi- 
lar to stellar systems in astrophysics. In Sees. 01 and El we 
shall discuss the kinetic theory of the HMF model and 
obtain the equivalent of the Vlasov and Landau equa- 
tions that are used to describe the dynamics of elliptical 
galaxies and globular clusters respectively. However, in or- 
der to facilitate the discussion and the comparison, it is 
useful to discuss first the dynamics of a one-dimensional 
barotropic fluid system with cosine interactions described 
by the Euler equations. In astrophysics, these equations 
describe the dynamics of barotropic stars. Stellar systems 
and barotropic stars are often treated in parallel due to 
their analogies (Binney & Tremaine 1987). In particular. 



it is possible to infer sufficient conditions of dynamical 
stability for spherical stellar systems from the dynamical 
stability of a barotropic star with the same density distri- 
bution. This constitutes the Antonov first law. Therefore, 
it is also of interest to develop this parallel in the case of 
the HMF model. To have a similar vocabulary, the sys- 
tems considered in this paper will also be called "stellar 
systems" and "gaseous stars" although they are only one- 
dimensional and correspond to a cosine interaction. 



3.1 Euler equations and energy functional 

We consider a gaseous system described by the Euler equa- 
tions 



9p d , . ^ 



du 
'dt 



de 



1 dp 
pde 



'89' 



(66) 



(67) 



where the potential is given by (jSJ. To close the equa- 
tions, we consider an arbitrary barotropic equation of state 
P = P{p)- We emphasise that these equations cannot be 
derived from the HMF model (Q) which rather leads to 
kinetic equations like the Vlasov equation. However, we 
shall see that there is a close connection between the sta- 
tionary states of the Vlasov and the Euler equations and 
that the limits of dynamical stability are the same in the 
two systems. Thus, the study of the Euler equation (which 
is simpler than the Vlasov equation) brings many infor- 
mation about the stability of stationary states of the HMF 
model with respect to the Vlasov equation even if the Eu- 
ler system does not describe dynamically the HMF model. 

It is straightforward to verify that the energy func- 
tional 



W = 



(68) 



is conserved by the Euler equations [W = 0). The first 
term is the internal energy, the second the potential en- 
ergy and the third the kinetic energy associated with the 
mean motion. The mass is also conserved. Therefore, a 
minimum of W at fixed mass determines a stationary so- 
lution of the Euler equations which is formally nonlinearly 
dynamically stable in the sense of Holm et al. (1985). We 
are led therefore to consider the minimization problem 

Min {W[p,u] I M[p] = M}. (69) 



3.2 First variations: the condition of hydrostatic 
equilibrium 

Cancelling the first order variations of Eq. H68() , we obtain 
u — and the condition of hydrostatic equilibrium 



dp 

de 



d$ 
d0' 



(70) 
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Therefore, extrema of W correspond to stationary solu- 
tions of the Euler equations (|66(l - (|67|l . On the other hand, 
combining the condition of hydrostatic equilibrium l|70|) 
and the equation of state p — p{p), we get 



(71) 



so that p is a function of that we note p = p(^). Using 
Eq. (0, we find that 

p = p{B cose), (72) 
where B is determined by the implicit equation 



B = / p{B cose') cos 9' de' 

27r Jo 



(73) 



Again, _B = is a solution of this equation characterizing 
a homogeneous phase p = M /iir. To determine the point 
of bifurcation to the inhomogeneous phase, we expand Eq. 
H73(l around _B = 0. Then, we find that cluster solutions 
appear when 



k dp , ^ 



(74) 



Using the condition of hydrostatic balance (|70|l , this can 
be rewritten 



(75) 



where Cg — {dp/dp)^^'^ is the velocity of sound in the ho- 
mogeneous phase where p = AI /2'k. 



3.3 Second variations: the condition of nonlinear 
dynamical stability 



phase is formally nonlinearly dynamically stable with re- 
spect to the Euler equations when 



kM 
1^' 



(78) 



and dynamically unstable otherwise. According to Eq. H75|l . 
the onset of dynamical instability coincides with the point 
where the clustered phase appears. 

The stability of the clustered phase can be investigated 
by solving the eigenvalue equation (|77|l for a specified 
equation of state p{p). This equation is the counterpart of 
equation (222) of Chavanis (2003a) for 3D self-gravitating 
gaseous spheres. 



3.4 The condition of linear stability: Jeans-like 
criterion 

We now study the linear dynamical stability of a station- 
ary solution of the Euler equation. We consider a small 
perturbation around a stationary solution of Eqs. (|66() - 
(|57|l and write p = p + Sp, u = Su etc... The linearized 
equations for the perturbation are 



dSp d 



+ -7TJ\{P6U) 0, 



dt de 



(79) 



d5u d d5(I> d<l> 

s^ie) = — / cos(6i - e')Sp{e')de'. (si) 

Jo 

Writing the time dependence in the form Sp ^ e^^,..., we 
get 



0, 



(82) 



The second variation of W due to variation of the velocity 
is trivially positive. The second variation of W due to 
variation of p is 

5^w = ^J 5p5me + J ^^{Sp)^de, (76) 

which must be positive for stability. Using the same pro- 
cedure as in Sec. 12.51 we find that the eigenvalue equation 
determining the stability of the solution is now 

(77) 

and that the condition of stability is A < (this yields a 
maximum of — W). For the uniform solution p = M/2'k, 
we can repeat exactly the same steps as in Sec. 12.61 since 
p'{p) is a constant which plays the role of T in the ther- 
modynamical analysis. Therefore, we find that the uniform 



d dS<P d<P 

Introducing the notation (|50|l . the continuity equation can 
be integrated into 

Xq + p5u ^ 0, (84) 

where we have imposed Su{0) — Su(2'k) = 0. Substitut- 
ing this relation in Eq. I|83|) and using the condition of 
hydrostatic equilibrium H70() . we finally obtain 

(85) 

This equation is the counterpart of the Eddington equa- 
tion of pulsations for a barotropic star (see also Eq. (224) 
of Chavanis 2003a). We note that Eqs. lO and ijSSJ) coin- 
cide for the neutral point A = 0. Therefore, the conditions 
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of linear stability and formal nonlinear dynamical stabil- 
ity coincide. The same conclusion holds for 3D barotropic 
stars (Chavanis 2003a). 

Considering the uniform phase p = M /2n and follow- 
ing a method similar to that developed in Sec. 12.61 we find 
that the most destabilizing mode (n = 1) is 



27rA 



5p = aiCOs9e , Su = ai sin 6*6 , 

M 



where the growth rate is given by 

47r ' 



(86) 



(87) 



When < kM/An, then A — ±vA^ and the perturbation 
grows exponentially rapidly (unstable case). When > 
kM/AiT, then A = ±^^/-\^ and the perturbation oscillates 
with a pulsation uj = V— A^ without attenuation (stable 
case). Therefore, the uniform phase is linearly (and also 
formally nonlinearly) dynamically stable with respect to 
the Euler equations when 



kM 
1^' 



(88) 



and linearly dynamically unstable otherwise. 

3.5 Particular examples 

3.5.1 Isothermal gas 

For an isothermal gas, we have 

P = PT, cl = T, 

and 



(89) 



W = T / pin 



p-d0, (90) 



where the temperature T is uniform. We note that the 
energy functional (|90|l of an isothermal gas coincides with 
the Boltzmann free energy Fb[p] — E[p] — TSb[p] of a 
A^-body system in the canonical ensemble, see Eqs. (|45ll 
and (|46|l of Sec. 12.51 This remark also holds for 3D self- 
gravitating systems (Chavanis 2002b). The pulsation equa- 
tion H85|l becomes 



d f 1 dq 



d0\pd0j 2ttT Jq 



qyu ) cos( 



(91) 



which can be connected to Eq. (|56|) . According to the cri- 
teria (|78(l - (|88|l . the uniform phase is formally nonlinearly 
dynamically stable for 



T>Tr. 



kM 
1^' 



(92) 



and linearly dynamically unstable otherwise. This crite- 
rion (or more generally the criterion 1)88(1 ) can be regarded 







B=0 - - 




/ 

/ 


UNSTABLE 


/ 

o 1 AbLb/ 

/ 

/ 

/ 


\ 


STABLE 


■ / 




B>0 


B=o 

1 

1 







1 2 3 4 5 

Fig. 9. Growth rate and pulsation period of an isothermal gas 
as a function of the temperature. 



as the counterpart of the Jeans instability criterion in 
astrophysics (Binney & Tremaine 1987). We emphasize, 
however, an important difference. In the case of 3D self- 
gravitating systems, the Jeans criterion selects a critical 
wavelength Aj (increasing with the temperature) above 
which the system is unstable against gravitational col- 
lapse. In the present context, where the interaction is trun- 
cated to one Fourier mode n = 1, the criterion ((92|l selects 
a critical temperature below which the system is unstable. 
The generalization of the Jeans instability criterion for an 
arbitrary binary potential of interaction in D dimensions 
is discussed in Appendix ICl and in Chavanis (2004b). This 
generalization clearly shows the connection between the 
HMF model and 3D self-gravitating systems. 

According to Eq. H87|) , the relation between A and the 
temperature T is 

kM 



47r 



(93) 



For T <T^, the growth rate is A = (T^ - T)^^^ and for 
T > Tc, the pulsation is w = (T - Tc)^/^. Following the 
preceding remark, we stress that, in the present context, 
A and to only depend on the temperature T, while for a 
3D gravitational gas, they depend on the wavelength of 
the perturbation (Binney & Tremaine 1987). Here, the 
unstable mode is fixed to n = 1. 

Considering now the clustered phase and using a per- 
turbative approach similar to that of Appendix^for T 
T~ (not detailed), we find that the pulsation is given by 
oj = \Jl{Tc — T). For smaller temperatures, the eigen- 
value equation can be solved numerically and the results 
are shown in Fig. 



3.5.2 Polytropic gas 

For a polytropic gas, we have 

p = Kp'^, 7 



(94) 



12 



P.H. Chavanis, J. Vatteville and F. Bouchet: Dynamics and thermodynamics of the HMF model 



where K is the polytropic constant and n is the polytropic 
index. For n +cxi, we recover the isothermal case with 
7=1 and K — T. For that reason K is sometimes called 
a polytropic temperature. The energy functional 1)68(1 can 
be written 



6 



p-de. (95) 



We have added a constant term (proportional to the total 
mass) in the polytropic energy functional so as to re- 
cover the isothermal energy functional ((90|1 for n — > +00. 
Under this form, we note that the energy functional of a 
polytropic gas has the same form as the Tsallis free en- 
ergy F^[p] = E[p] — KS^[p] where 7 plays the role of the 
g-parameter and K the role of a generalized temperature. 
The same remark holds for 3D self-gravitating systems. 
However, this resemblance is essentially fortuitous and the 
mark of a thermodynamical analogy (Chavanis 2003a; Cha- 
vanis & Sire 2004b). 

If we define the local temperature by p{9) = p{d)T{d), 
we obtain T{9) = Kpififl'' and cliO) = -fT{e). We note 
that, for a polytropic distribution, the kinetic temperature 
T{6) usually depends on the position while the polytropic 
temperature K is uniform as in an isothermal gas. How- 
ever, in the uniform phase T — Kp^/^ is a constant that 
can be called the temperature of the polytropic gas. The 
velocity of sound in the homogeneous phase is 



1 



= ,T^K,p^-^^Kt^(^ 
n \2iT 



l/n 



(96) 



The condition of dynamical stability ((78|l - 188|l can be writ- 
ten 



K>Kn 



kM 



47r 1 



or, equivalently. 



1 



kM 



l/r. 



(97) 



(98) 



For 7 > 1 (i.e., n > 0), the critical temperature T„ is 
smaller than the corresponding one for an isothermal gas 
Tc = Too, i-e. the instability is delayed. For 7 < 1 (i.e., 
n < 0), the instability is advanced. Similar results are 
obtained for 3D gravitational systems (Chavanis & Sire 
2004b). According to Eq. l(57jl . the relation between A and 
the kinetic temperature T is 



A2 = 7(T„-r). 



3.6 The local Euler equation 



(99) 



We can consider a simplified problem where the potential 
in Eq. H67|l is fixed to its equilibrium value <P = BcosO. 
In that case, we get the local Euler equation 



du 
'dt 



du 



Idp 

'~pde 



Bsine. 



(100) 
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Fig. 10. Characteristics of the forced Burgers equation 11U2II 
showing the appearance of shocks and chevrons. 



The stationary solutions are given by H72|) . The linear 
stability of a stationary solution amounts to solving the 
Sturm-Liouville problem 



d / p'{p) dq 
d0\ p de 



A2 



P 



-q- 



(101) 



This problem will be considered in Sec. l6.3l for the isother- 
mal equation of state. 

In the pressureless case, we get the forced Burgers 
equation 



du du 



(102) 



It can be solved by the method of characteristics, writing 
the equation of motion of a particle as 



dt^ 



B sin e. 



(103) 



The trajectories 6{t) can then be expressed in terms of el- 
liptic functions. The dynamics of the forced Burgers equa- 
tion is interesting as it develops "shocks" and "chevrons" 
(caustics) singularities (see Fig. [Tn|) . The Burgers equa- 
tion also appears in cosmology to describe the formation 
of large-scale structures in the universe (Vergassola et al. 
1994). A detailed description of this dynamics is given by 
Barre et al. (2002) in the context of the repulsive HMF 
model. In that case, the forced Burgers equation H102|l 
with sin(26') instead of sinS models the short time dy- 
namics of the Hamiltonian A^-body system. For the at- 
tractive HMF model the short time dynamics can be 
modelled by the non-local Euler equation 1(67(1 with zero 
pressure p — 0. This is solution of the Vlasov equation (see 
Sec.lU with f{e,v,t) = p{0,t)S{v - u{0,t)). This single- 
speed solution is valid until the first shock. However, the 
connection with the local Euler equation ((102(1 is not clear 
in that case because the homogeneous phase is unstable 
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which precludes the possibility of deriving p()2|l from l|67|l 
as is done in Barre et al. (2002) in the repulsive case. Note 
also that in the attractive HMF model (ferromagnetic) we 
just have one cluser while the repulsive HMF model (anti- 
ferromagnetic) shows a bicluster. 



4 Violent relaxation, metaequilibrium states 
and dynamical stability of collisionless stellar 
systems 

We now come back to the HMF model defined by the 
Hamilton equations and develop a kinetic theory by 
analogy with stellar systems. In particular, we emphasize 
the importance of the Vlasov equation and the concept of 
violent relaxation introduced by Lynden-Bell (1967). 

4.1 Vlasov equation and H-functions 

For systems with long-range interactions, the relaxation 
time toward the statistical equilibrium state (|18|l is larger 
than Ntu, where to is the dynamical time (Chavanis 
2004b). Accordingly, for iV — > +oo, the relaxation time is 
extremely long and, for timescales of physical interest, the 
evolution of the system is essentially collisionless. More 
precisely, for t <C Ueiax and N +oo (Vlasov limit), the 
time dependence of the distribution function is governed 
by the Vlasov equation 



dt de de dv 



(104) 



which has to be solved in conjunction with Eq. ©. This 
system of equations is similar to the Vlasov-Poisson sys- 
tem describing the dynamics of elliptical galaxies and other 
collisionless stellar systems in astrophysics. Starting from 
an unstable initial condition, the HMF model ^ will 
achieve a metaequilibrium state (on a coarse-grained scale) 
as a result of phase mixing and violent relaxation (Binney 
& Tremaine 1987). This metaequilibrium state is a par- 
ticular stationary solution of the Vlasov equation. Since 
it results from a complex mixing, it is highly robust and 
nonlinearly dynamically stable with respect to collision- 
less perturbations. The process of violent relaxation and 
the convergence of the distribution function toward a sta- 
tionary solution of the Vlasov equation has been illus- 
trated numerically by Yamaguchi et al. (2004) for the 
HMF model. This is similar to the violent relaxation of 
stellar systems in astrophysics and 2D vortices in hydro- 
dynamics (see Chavanis 2002a). 

One question of great importance is whether we can 
predict the metaequilibrium state achieved by the system 
as a result of violent relaxation. Lynden-Bell (1967) has 
tried to make such a prediction by resorting to a new 
type of statistical mechanics accounting for the conser- 
vation of all the Casimirs imposed by the Vlasov equa- 
tion. This theory was developed for the gravitational in- 
teraction, but the general ideas and formalism apply to 
any system with long-range interactions described by the 



Vlasov equation. In the non-degenerate limit, he predicts 
a Boltzmann distribution of the form / ~ e~^'^ where the 
individual mass of the particles does not appear. Lynden- 
Bell (1967) also understood that his statistical prediction 
is limited by the concept of incomplete relaxation. The 
system tries to reach the most mixed state but, as the 
fluctuations become weaker and weaker as we approach 
equilibrium, it can settle on a stable stationary solution 
of the Vlasov equation which is not the most mixed state. 
In order to quantify the importance of mixing, Tremaine 
et al. (1986) have introduced the concept of H-functions 



C{f)d9dv, 



(105) 



where C is an arbitrary convex function. The H-functions 
c_alculated with the coarse-grained distribution function 
/ increase as a result of phase mixing in the sense that 
S[J{0,v,t)] > 5'[7(6',w,0)] for t > where it is assumed 
that, initially, the system is not mixed so that f{0, v, 0) = 
f{9,v,0). This is similar to the i/-theorem in kinetic the- 
ory. However, contrary to the Boltzmann equation, the 
Vlasov equation does not single out a unique functional 
(the above inequality is true for all i?-functions) and the 
time evolution of the i/-functions is not necessarily mono- 
tonic (nothing is implied concerning the relative values of 
H{t) and H{t') for t,t' > 0). On the other hand, any 
stationary solution of the Vlasov equation of the form 
f ^ f(e) with /'(e) < extremizes a if- function at fixed 
mass and energy. If, in addition, it maximizes S at fixed E, 
M, then it is nonlinearly dynamically stable with respect 
to the Vlasov equation. In astrophysics, such distribution 
functions depending only on the energy describe spherical 
stellar systems. This is a particular case of the Jeans the- 
orem (Binney & Tremaine 1987). For a ID system such as 
the HMF model, this is the general form of inhomogeneous 
stationary solutions of the Vlasov equation. Therefore, we 
expect that the H-functions will increase during violent 
relaxation until one of them (non-universal) reaches its 
maximum value at fixed mass and energy when a sta- 
tionary solution of the Vlasov equation is reached (this is 
not necessarily the case in astrophysics since the system 
can reach a steady state that does not depend only on 
energy). Note that the Boltzmann and the Tsallis func- 
tional are particular H-functions (not thermodynamical 
entropies in that context) associated with particular sta- 
tionary solutions of the Vlasov equation: isothermal stellar 
systems and stellar polytropes (Chavanis & Sire 2004b). 
All these ideas, first developed for stellar systems, apply 
to other systems with long-range interactions such as the 
HMF model. 



4.2 Nonlinear dynamical stability criterion for the 
Vlasov equation 

The theory of violent relaxation explains how a collision- 
less system out of mechanical equilibrium can reach a 
steady solution of the Vlasov equation on a very short 
timescale due to long-range interactions and chaotic mix- 
ing. Since this metaequilibrium state is stable with respect 



14 



P.H. Chavanis, J. Vatteville and F. Bouchet: Dynamics and thermodynamics of the HMF model 



to coUisionless perturbations, it is of interest to determine 
a criterion of formal nonlinear dynamical stability for the 
Vlasov equation. For the HMF model, this has been con- 
sidered by Yamaguchi et al. (2004) using the Casimir- 
Energy method. We shall propose another derivation of 
the stability criterion which uses a formal analogy with the 
thermodynamical analysis developed in Sec. |21 and which 
is also applicable to the clustered phase. This approach is 
similar to the one developed by Chavanis (2003a) for 3D 
stellar systems. 

Let us introduce the functional S = — J C{f)d9dv 
where C (/) is a convex function. This functional is a par- 
ticular Casimir so it is conserved by the Vlasov equation. 
The energy E and the mass M are also conserved. There- 
fore, a maximum of S at fixed mass M and energy E deter- 
mines a stationary solution f{9, v) of the Vlasov equation 
that is nonlinearly dynamically stable. We are led there- 
fore to consider the maximization problem 

Max {S[f] I E[f] = E, M[f] - A/}. (106) 

We also note that F[f] = E[f] - TS[f] (where T is a 
positive constant) is conserved by the Vlasov equation. 
Therefore, a minimum of F at fixed mass M is nonlinearly 
dynamically stable with respect to the Vlasov equation {F 
is called an energy-Casimir functional). This corresponds 
to the formal stability criterion of Holm et al. (1985). This 
criterion can be written 

Min{i^[/] \M[f]=M}. (107) 

To study the nonlinear dynamical stability of collision- 
less stellar systems, we are thus led to consider the max- 
imization problems 1)106(1 and 1107|l . These are similar to 
the conditions of thermodynamical stability (|15|l and H16|l 
but they involve a more general functional S[f] than the 
Boltzmann entropy. In addition, they have a completely 
different interpretation since they determine the nonlin- 
ear dynamical stability of a steady solution of the Vlasov 
equation, not the thermodynamical stability of the statis- 
tical equilibrium state. Due to this formal resemblance, we 
can develop a thermodynamical analogy (Chavanis 2003c) 
and use an effective thermodynamical vocabulary to inves- 
tigate the nonlinear dynamical stability of a coUisionless 
stellar system. In this analogy, S plays the role of an effec- 
tive entropy, T plays the role of an effective temperature 
and F plays the role of an effective free energy. The cri- 
terion H106|l is similar to a condition of microcanonical 
stability and the criterion (|107|l is similar to a condition 
of canonical stability. 

We also note that the stability criterion (|106|l is con- 
sistent with the phenomenology of violent relaxation. In- 
deed, the H-functions increase on the coarse-grained scale 
while the mass and the energy are approximately con- 
served. Therefore, the metaequilibrium state is expected 
to maximize a certain H-function (non-universal) at fixed 
mass and energy. In that interpretation, / has to be viewed 
as the coarse-grained distribution function /, not the dis- 
tribution function itself. The point is that during mixing 
Df/Dt and the 7J-functions S[f] increase. Once it 



has mixed Df/Dt = so that S[f] = 0. Since f{e,v,t) 
has been brought to a maximum /q(^,w) of a certain H- 
function and since S[f] is conserved (after mixing), then 
/q is a nonlinearly dynamically stable steady state of the 
Vlasov equation. 

4.3 First variations: stationary solutions of the Vlasov 
equation 

Introducing Lagrange multipliers as in Sec. 12.31 the criti- 
cal points of the variational problems (|106|l and H107|) are 
given by 

C7'(/) = -/3e-a, (108) 

where e = ^ -f ^ is the energy of a particle. Since C is 
a monotonically increasing function of /, we can inverse 
this relation to obtain 

f = F{Pe + a), (109) 

where F{x) = (C")~^(— x). We can check that any DF 
/ — /(e) is a stationary solution of the Vlasov equation 
(|104() . From the identity 

/'(6) = -/3/C"(/), (110) 

resulting from Eq. H108() . we see that /(e) is a monotonic 
function of the energy. Assuming that /(e) is decreasing, 
which is the physical situation, imposes /3 = 1/T > 0. 

We note also that for each stellar system with / = 
/(e), there exists a corresponding barotropic star with 
the same equilibrium density distribution. Indeed, defin- 
ing the density and the pressure hy p = fdv = p{^) , 
p — J/^^ fv'^dv = p{^), and eliminating the potential 
^ between these two expressions, we find that p = p{p). 
Writing explicitly the density and the pressure in the form 

P = 2/^"°°m + a)^^^-L^de, (111) 

p + CO 

P = 2 F{Pe + a)[2{e--^)]'/^de, (112) 

J* 

and taking the 6'-derivative of Eq. H112I) . we obtain the 
condition of hydrostatic equilibrium l(7n|l . 

Due to the analogy between stellar systems and barotropic 
stars, it becomes possible to use the results obtained in 
Sec.Olto study the stationary solutions of the Vlasov equa- 
tion (|104(l . In particular, the transition from homogeneous 
{B = 0) to inhomogeneous {B ^ 0) solutions is again 
given by the criterion (|75|l . Now, in the case of stellar sys- 
tems, it is more relevant to express this criterion in terms 
of the distribution function. Using the identity ((299|l . the 
criterion (|75|l determining the appearance of the clustered 
phase is equivalent to 

i+^r/>id.<o. (113) 

2 J-oo V 

We will soon see how this quantity is related to the dielec- 
tric function of a gravitational plasma. 
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4.4 Second variations: the condition of nonlinear 
dynamical stability 

We shall investigate the formal nonlinear dynamical sta- 
bility of stationary solutions of the Vlasov equation by us- 
ing the criterion (|107|l . This criterion is less refined than 
the criterion H1U()|) because all solutions of 1)1071) are so- 
lution of (|106|l . but the reciprocal is wrong in general. In 
particular, for long-range interactions, the optimization 
problems H1U6|) and H1U7|) may not coincide. In thermo- 
dynamics, this corresponds to a situation of ensemble in- 
equivalence (see Bouchet & Barre 2004). Therefore, the 
criterion l|107l) can only give a sufficient condition of non- 
linear dynamical stability for stationary solutions of the 
Vlasov equation of the form / = /(e) with /'(e) < 0. This 
corresponds to the criterion of formal nonlinear dynami- 
cal stability given by Holm et al. (1985). The more refined 
criterion (|106|l has been introduced by Ellis et al. (2002) 
in 2D hydrodynamics (for the 2D Euler-Poisson system) 
and applied by Chavanis (2003a) in stellar dynamics (for 
the Vlasov-Poisson system). 

To obtain a managable criterion of dynamical stabil- 
ity, we use the same procedure as the one developed in 
Chavanis (2003a). We shall not repeat the steps that are 
identical. We first minimize the functional F[f] at fixed 
temperature and density p{9). This gives an optimal dis- 
tribution f*{e,v), determined by C"(/*) = -(3^ - X{9), 
which depends on the density p{9) through the Lagrange 
multiplier X{0). Then, after some manipulations, we can 
show that the functional F[p] — F[/*] can be put in the 
form 



F = 



p<Pde 







(114) 



where p{p) is the equation of state determined by C(/) 
according to Eqs. Ullll) and (|112|l . We now need to min- 
imize F[p] at fixed mass. To that purpose, we just have 
to observe that F[p] corresponds to the energy functional 
H68I) of a barotropic gas with u = 0. Therefore, the can- 
cellation of the first variations of Eq. pi4|l returns the 
condition of hydrostatic equilibrium (|70|l and the positiv- 
ity of the second variations leads to the stability criterion 
A < linked to the eigenvalue equation H77|l . Therefore, 
the criterion of formal nonlinear dynamical stability (|107|l 
for the Vlasov equation (stellar systems) is equivalent to 
the criterion of formal nonlinear dynamical stability l|69(l 
for the Euler equations (gaseous barotropic stars). 

Using the results of Sec. 01 we conclude that the uni- 
form phase is formally nonlinearly dynamically stable with 
respect to the Vlasov equation when 



kM 



(115) 



When the inequality is reversed, the uniform phase is a 
saddle point of F at fixed mass M and we shall see that it 
is linearly dynamically unstable. Using the identity 1)2991) . 
the nonlinear criterion H115() can be rewritten as 



+ 00 



fiv) 



dv > 0, 



(116) 



which was found by Yamaguchi et al. (2004) using a dif- 
ferent method. An advantage of the present approach is 
that this approach is also applicable to an inhomogeneous 
system. Indeed, the stability of the clustered phase can be 
investigated by solving the eigenvalue equation H77|) for 
the equation of state specified by the function C(/), and 
investigating the sign of A. 



4.5 About the Antonov first law 

As discussed previously, the criterion (|107|l providing a 
condition of nonlinear dynamical stability for a stellar sys- 
tem with / = /(e) and /'(e) < is equivalent to the crite- 
rion H69I) determining the nonlinear dynamical stability of 
a barotropic star with the same equilibrium density dis- 
tribution. On the other hand, we have already indicated 
that the criterion H107() is less refined than the criterion 
(|106)) which is believed to be the strongest criterion of 
nonlinear dynamical stability for stationary solutions of 
the Vlasov equation of the form f — f{e) with /'(e) < 0. 
In general, the criterion l|107l) just provides a sufficient 
condition of nonlinear dynamical stability. Thus, we can 
"miss" stable solutions if we use just (|l()7)l instead of 1)106)1 
[said differently, the set of solutions of ()107)) is included 
in 1)106)) ]. From these remarks, we conclude that "a stellar 
system is stable whenever the corresponding barotropic 
gas is stable" but the converse is wrong in general. This 
is the so-called Antonov first law in astrophysics (Binney 
& Tremaine 1987). Our approach provides an extension of 
the Antonov first law to the case of nonlinear dynamical 
stability (while the usual Antonov first law corresponds to 
linear dynamical stability). Furthermore, by developing a 
thermodynamical analogy, we have provided an original 
interpretation of the nonlinear Antonov first law in terms 
of "ensembles inequivalence" for systems with long-range 
interactions (Chavanis 2003a) . 

For 3D self-gravitating systems, we know that the en- 
sembles are not equivalent so that the criterion ()107)) is 
more restrictive than the criterion 1)106)) . In that case, 
()107)) is just a sufficient condition of nonlinear dynamical 
stability (definitely). For example, using the criterion 1)69)) . 
it can be shown that polytropic stars with index n < 3 are 
nonlinearly dynamically stable with respect to the Euler 
equations while polytropic stars with index 3 < n < 5 
are dynamically unstable (polytropes with index n > 5 
have infinite mass). Therefore, using the "canonical" cri- 
terion ()107)l . we can only deduce that stellar polytropes 
with index ti < 3 are nonlinearly dynamically stable with 
respect to the Vlasov equation. However, using the "mi- 
crocanonical" criterion ()106)l , we can prove that all stellar 
polytropes with index n < 5 are nonlinearly dynamically 
stable with respect to the Vlasov equation. Polytropes 
with index 3 < n < 5 lie in a region of "ensemble in- 
equivalence" (in the thermodynamical analogy) where the 
"specific heat" is negative (Chavanis 2003a). 

For the HA4F model considered in this paper, we be- 
lieve that the criteria ()106)l and 1)107)1 determine the same 
set of solutions so that the ensembles are equivalent in that 
case (no solution is forgotten by 1)107)) '). We have shown 
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in Sec. El that this is at least the case for isothermal dis- 
tributions. If wc take for granted that this equivalence 
extends to any functional (|105|l . we conclude that, for the 
HMF model, "a stellar system is stable if, and only if, 
the corresponding barotropic gas is stable". This would 
be the HMF version of the nonlinear Antonov first law. 
In that case, the stability limits obtained for the Euler 
equation in Sec.|3|can be directly applied to the Vlasov 
equation (they coincide). Clearly, it would be of great in- 
terest to derive general criteria telling when the ensembles 
are equivalent or inequivalent, depending on the form of 
the functional H1Q5(I and on the form of the potential of 
interaction u{r — r'). 



extremizes a H-function at fixed mass and energy, then 
it has a single maximum at v = 0. Therefore, ujr = 
according to Eq. H120I) and the criterion of linear dynam- 
ical stability H121|l coincides with the criterion of formal 
nonlinear dynamical stability (|116|l . 



4.7 Particular examples 

We shall now present explicit results for particular sta- 
tionary solutions of the Vlasov equation. We consider the 
case of isothermal stellar systems and stellar polytropes. 



4.6 The condition of linear stability 

We now study the linear dynamical stability of a spatially 
homogeneous stationary solution of the Vlasov equation 
described by / = /(«). This problem was first investigated 
by Inagaki & Konishi (1993) and Pichon (1994) and more 
recently by Choi & Choi (2003). We shall complete here 
their study. Writing the perturbation in the form 6f ^ 
^i(ne-ujt) g^j^j using standard methods of plasma physics, 
we obtain the dispersion relation 

t{n,u) = 1 + J(<5„,i - 5„,_i) / -l^dv = 0, 
2 J nv — uj 

(117) 

where e(n,w) is the dielectric function and the integral 
must be performed by using the Landau contour. For the 
destabilizing mode n — 1 {n — —I gives the same result), 
Eq. (|117|l reduces to 

1 + ^ / I^dv = 0. (118) 

2 J V — U! 

The condition of marginal stability is = where uji 
is the imaginary part of w = Wr + iuji. In that case, the 
integral in Eq. H118|) can be written as 

l+'^V [ -l^dv + in^f'iur) = 0, (119) 

2 J V — LOr 2 

where V denotes the principal value. Identifying real and 
imaginary parts, if follows that 

k f f'(v) , 

2 J V — OJr 
f'{iUr) = 0. (120) 

The second relation fixes the frequency of the perturbation 
and the first equation determines the point of marginal 
stability in the series of equilibria. The system is linearly 
dynamically stable if 

1 + f -^^dv > 0, (121) 

2 J V — UJr 

and linearly dynamically unstable otherwise. Note that 
f{v) does not need to be symmetrical. However, if f(v) 



4.7.1 Isothermal stellar systems 

We consider the H-function 

S = - J finfdOdv, (122) 

which is similar to the Boltzmann entropy H14() in thermo- 
dynamics. However, as explained in Sec. 14.11 its physical 
interpretation is different. Its maximization at fixed mass 
and energy determines a formally nonlinearly dynamically 
stable stationary solution of the Vlasov equation corre- 
sponding to the isothermal distribution function 

/ = Ae^f^'. (123) 

This distribution function has the same form (but a dif- 
ferent interpretation) as the statistical equilibrium state 
(fl^ of the A^-body system. 

The barotropic gas corresponding to the isothermal 
distribution function H123(l is the isothermal gas with an 
equation of state p — pT where T = 1/(3. Therefore, the 
velocity dispersion of an isothermal stellar system is 
equal to the velocity of sound = T in the corresponding 
isothermal gas. The functional H114|) takes the form 

F[p] J pme + T J pin pd9, (124) 

and the density is related to the potential according to the 
formula 

p = A'e-'^*, (125) 

which can be obtained by extremizing F[p] at fixed mass. 
We can also express the distribution function in terms of 
the density according to 

/ a \ 1/2 2 

/=(^£j p{0)e-^^. (126) 

According to what has been said in Sec. 14. 5! about the cor- 
respondance between stellar systems and barotropic stars, 
we conclude that the uniform phase of an isothermal stel- 
lar system (|123ll is formally nonlinearly dynamically stable 
with respect to the Vlasov equation if T > Tc and linearly 
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dynamically unstable otherwise. In terms of the dimen- 
sionless parameters 77 = and e 
of dynamical stability can be written 

f? < 1, e > 1- 



-j^f^, the conditions 



(127) 



They coincide with the conditions of thermodynamical 
stability (see Sec. Ill- 

For the Maxwellian distribution function (|126|l with 
uniform density, the dielectric function can be written 



where 



W{z) 



/2Tr 



-|-oo 



2 dx, 



(128) 



(129) 



is the M^- function of plasma physics (Ichimaru 1973). This 
is an analytic function in the upper plane of the complex 
z plane which is continued analytically into the lower half 
plane. Explicitly, 



W{z) ze 



dye' 



ITT 

i\l —ze 



(130) 



We look for solutions of the dispersion relation e(l, lo) — Q 
in the form uj — iX where A is real. First, we note that 



e(l,iA) = 1 



where we have defined the G-function 
G{x) ^ 



1 — -JiTxe^^ eilc{x) 



(131) 



(132) 



For a; — > 0, G{x) = 1 + ^/^Tx + .... For x - 
2x^(1 + A + ■•■)• For X -00, G{x) 



-00, G{x) 

2^x 



Therefore, the relation between A and T is given by 



(133) 



The case of neutral stability uj — Q corresponds to T = 
Tc (or ?7 = 1). The case of instabihty (A > 0) corresponds 
to T < Tc- The perturbation grows exponentially rapidly 
as 6f ^ e^^. The growth rate A is given by Eq. H133|) which 
can be explicitly written 



1- 

T 



— -Ae^Terfc 
2T 



For T ^ , we have 
A- 

and for T — > 0, we have 

A- 



kM 



47r I Tr. 



1/2 



(134) 



(135) 



(136) 



0.5 



-0.5 



X>0 


stellar 
system 







0.5 



1.5 

T/T, 



2.5 



Fig. 11. Growth rate and decay rate of isothermal stellar sys- 
tems and isothermal stars as a function of the temperature in 
the framework of the HMF model. 



The first term in Eq. lfK!H|l can be deduced directly from 
Eq. H118|l by using the distribution function f{v) ^ p S{v) 
valid at T = and integrating by parts. The case of 
stability (A < 0) corresponds to T > Tc- The perturba- 
tion is damped exponentially rapidly as Sf ^ e^''* where 
7 — —A. This is similar to the Landau damping in plasma 
physics, except that here there is no pulsation (a;,- = 0). 
By contrast, in plasma physics, the pulsation tOr is much 
larger than the damping rate 7. The damping rate 7 = — A 
is given by 



V^F( ^/-7 



where we have defined the F-function 

1 



F{x) 



1 + V^oi 



rfc(— x) 



(137) 



(138) 



such that F{x) ~ G{-~x). For x 
For X -00, F{x) ^ 2x^(1 -f 

^(^)~ 2;be-'. Explicitly, 



0, F{x) = 1-^x- 
+ ...). For x - 



-00, 



— 7e2T'erfc 
2T ' 



For T^T+, we have 

and for T +00, we have 

7- V2TlnT. 



= 0. (139) 



(140) 



(141) 



Obviously, the relation H134|l between the growth rate 
and the temperature of an isothermal stellar system is 
different from the corresponding relation H93|l valid for an 
isothermal gas (they coincide only at T = 0). A similar 
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distinction is noted in the case of 3D self-gravitating sys- 
tems. In particular, Fig. 1111 can be compared with Fig. 
5.1 of Binney & Tremaine (1987). Moreover, in the un- 
stable regime, a gaseous medium supports sound waves 
with pulsation uj = {T — Tc^^^ that are not attenuated. 
By contrast, in a stellar medium at T < Tc, there exists 
solutions with no wave {tOr = 0) for which the perturba- 
tion is damped exponentially. Other solutions with 7^ 
probably exist but they are more difficult to investigate 
analytically. This is left for a future study. 

The relations (|134|l and H139|) have been obtained pre- 
viously by Choi & Choi (2003) using a sHghtly different 
approach. Our derivation emphasizes the close link with 
results in plasma physics. In addition, our formalism will 
be used in Sec. 15.41 to show that the damping rate 7 of 
a perturbation is equal to the exponential decay of the 
time auto-correlation function of the force. The general- 
ization of these results for an arbitrary form of long-range 
potential is given in Chavanis (2004b). 



nonlinearly dynamically stable stationary solutions of the 
Vlasov equation called stellar polytropes. The fact that 
the criteria IjlOGI) and (|107|1 of nonlinear dynamical stabil- 
ity are similar to criteria of generalized thermodynamical 
stability is the mark of a thermodynamical analogy (Cha- 
vanis 2003c, Chavanis & Sire 2004b). 

Stellar polytropes are described by the distribution 
function 



/ 



(143) 



obtained from Eq. I|1U8|) . When the term in brackets is 
negative, the distribution function is set equal to / = 0. 
The index n of the polytrope in one dimension is related to 
the parameter q by the relation (Chavanis & Sire 2004a) 



1 



(144) 



4.7.2 Stellar polytropes 

We consider the H function 

Sg = --^J{P-f)dedv, (142) 

where g is a real number. This functional has been intro- 
duced by Tsallis (1988) in non-extensive thermodynam- 
ics. The aim was to develop a generalized thermodynam- 
ical formalism to describe quasi-equilibrium structures in 
complex media that are not described by the Boltzmann 
distribution. In this sense, Sq[f] is interpreted as a gen- 
eralized entropy and its maximization as a condition of 
(generalized) thermodynamical stability. In the context of 
Vlasov systems, we shall rather intcrprete Sq[f] as a par- 
ticular if- function (see Sec. l4.1|l and its maximization as a 
condition of nonlinear dynamical stability ^ . Its maximiza- 
tion at fixed mass and energy leads to a particular class of 

^ If we were to apply Tsallis thermodynamics in the con- 
text of violent relaxation, we would need to introduce the fine- 
grained distribution of phase levels p{0, v, rj) [which is the rele- 
vant probability field in that context] and replace the Lynden- 
Bell entropy Sl.b.[p] = — J pin pdddvd-q (Lynden-Bell 1967) 

by Sq[p] — ^ /(p' ~ p)d9dvdrj as suggested in Brands et 

al. (1999). In that case, Sq[p\ can be regarded as a generalized 
entropy trying to take into account non-ergodicity and lack of 
complete mixing in coUisionless systems with long-range inter- 
actions. In that point of view, the parameter q measures the 
efficiency of mixing (g = 1 if the system mixes well which is 
implicitly assumed in Lynden-Bell's statistical theory). In the 
two-levels approximation / G {0, 1} and in the dilute limit 
/ <C 1, Sq[p\ can be expressed as a functional of the coarse- 
grained distribution function f = J prjdri = p of the form 

SqlJ] = --^JiT- J)dedv. In this particular limit, Sq[f] 
can be interpreted as a thermodynamical entropy generaliz- 
ing S[f] = — J f In fdOdv which is a particular case of the 
Lynden-Bell entropy for two levels in the dilute limit (see Cha- 
vanis 2004c). In conclusion, Tsallis functional Sq[p] expressed 
in terms of p(6,v,rj) is an entropy but Tsallis functional Sq[f] 



For n — !■ -|-oo (or q —f 1), we recover the isothermal dis- 
tribution function H123|l and the H- function HI22(I . There- 
fore, Tsallis functional connects continuously isothermal 
and polytropic distributions. Physical polytropic distri- 
bution functions (see Chavanis & Sire 2004a) have (3 > 
and 9 > 1 (i.e. n > 1/2) or 1/3 < g < 1 (i.e. n < -1). 

The barotropic gas corresponding to the polytropic 
distribution function l|143|) is the polytropic gas 

p = Kp^, 7=1 + 1, (145) 



expressed in terms of f{9, v) is either a _ff-function (dynamics) 
or the reduced form of entropy Sq[p\ (thermodynamics) for 
two levels in the dilute limit. In any case, it is not clear why 
non-ergodic effects could be encapsulated in the simple func- 
tional Sq[p] introduced by Tsallis. Tsallis entropy is "natural" 
because it has mathematical properties very close to those pos- 
sessed by the Boltzmann entropy and it is probably relevant to 
describe a certain type of mixing and non-ergodic behaviour 
as in the case of porous media and weak chaos (it may be seen 
as an entropy on a fractal space-time). However, many other 
types of non-ergodic behaviour can occur and other functionals 
S = — J C{p)d6dvdri could be relevant. Observation of stel- 
lar systems, 2D vortices and quasi-equilibrium states of the 
HMF model resulting from incomplete violent relaxation do 
not favour Tsallis distributions in a universal manner. Other 
distributions can emerge. In fact, we must give up the hope to 
predict the metaequilibrium state in case of incomplete relax- 
ation. We must rather try to construct stable stationary so- 
lutions of the Vlasov equation in order to reproduce observed 
phenomena. The H-functions (11051 can be useful in that con- 
text. An alternative procedure is to keep the Lynden-Bell form 
of entropy unchanged and develop a dynamical theory of vio- 
lent relaxation (Chavanis et al. 1996, Chavanis 1998) in order 
to take into account incomplete mixing through a variable dif- 
fusion coefficient related to the strength of the fluctuations. In 
that case, non-ergodicity is explained by the decay of the fluc- 
tuations of $ driving the relaxation, not by a complex structure 
of phase-space. Generalized entropies are not necessary in that 
case. 



P.H. Chavanis, J. Vatteville and F. Bouchet: Dynamics and thermodynamics of the HMF model 



19 



with the polytropic constant 



K 



K =-- 



V2A 



r{i/2)r{n- 



r{n + i) 



1/2)| 



- 1 /n 



n > 



(146) 



1 



V2A 



r{i/2)ri~n) 

r{l/2-n) 



-1/n 



n<~l (147) 



where A = {\q — l\P/qy^^'^~^\ In the present context, 
the polytropic constant K is related to the Lagrange mul- 
tiplier /3. Therefore, K and Tq = (3~^ play the role of 
effective temperatures (see Chavanis & Sire 2004b for a 
more detailed discussion). For a polytropic distribution, 
the functional (|114|l takes the form 



F[p] P'^dO 



K 



(p^ ~ p)de, (148) 



and the relation between the density and the potential is 



A 



(149) 



Comparing Eqs. H143|l and H149|l . we note that a polytropic 
distribution with index q in phase space yields a polytropic 
distribution with index 7 = l + 2[q—l) / {q-\-\) in physical 
space. In this sense, polytropic laws are stable laws since 
they keep the same structure as we pass from phase space 
/ = /(e) to physical space p = p{'P) as noticed in Chava- 
nis (2004d). This is probably the only distribution enjoy- 
ing this property. Similarly, comparing H142() -(|7I) and (|148|l 
the "free energy" in phase space F[f] = E[f] - ToSqlf] 
(where Tq = 1/(3) becomes F[p\ — E[p\ -KS-y[p\ in physi- 
cal space. Morphologically, the polytropic temperature K 
plays the same role in physical space as the temperature 
Tq = 1/P in phase space. 

We can express the distribution function in terms of 
the density according to 



/4 



l/n _ 



vy2 



{n+l)K 



n-1/2 



(150) 



with 



V2. 



r(l/2)r(n+ 1/2) 
r{n+l) 



[X(r^+l)]l/^ n>- (151) 



^ = V2^^i(ffl-^[-X(n + l)]V^ „<-i. (152) 

Introducing the kinetic temperature (velocity dispersion) 
T{e) = (w^) = p{9)/p{9) = i4:p(6')i/", this can be rewrit- 
ten 
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1 - 



1-2/2 



(n + l)T(0) 



ri-1/2 



with 



Bn — 



r(n + 1) 



r(n+ l/2)(n+ 1)1/2^ 
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(153) 
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£=87iE/kM 

Fig. 12. Bifurcation diagram of stellar polytropes that are 
stationary solutions of the Vlasov equation. The homogeneous 
phase is nonlinearly dynamically stable for e > Scrit = 1/7 
where 7 = 1 + l/n and n = i H — It becomes linearly dy- 
namically unstable for e < ecrit w?iere the branch of clustered 
states (represented schematically by a line) appears. Isother- 
mal stellar systems correspond to (g = 1, n = 00, 7 = 1). In 
ordinate, T is defined by K{M /2'kY^" . In the homogeneous 
phase, it represents the kinetic temperature of a polytropic 
stellar system. Note that is a monotonic function of the La- 
grange multiplier /3 so that the curve can be viewed as a series 
of equilibria of polytropic distributions. 



Bn — 



r{i/2- 



r(-n)[-(rH-l)]i/2 



n < -1. 



(155) 



Equation p53|l is the counterpart of Eq. H126(l for isother- 
mal systems. For n > 1/2, the distribution / = for 

\v\ > Vmax = \/2{n + 1)T. 

According to what has been said in Sec. 14. 41 about the 
correspondance between stellar systems and barotropic 
stars, we conclude that the uniform phase of a polytropic 
stellar system ()143ll with index n is formally nonlinearly 
dynamically stable with respect to the Vlasov equation if 
K > Kn or T > T„ and linearly dynamically unstable 
otherwise. It can be useful to introduce the dimensionless 
parameter 77 = kM/AnT where T = Kp^/^ is the kinetic 
temperature in the homogeneous phase where p = M /2tt. 
For n +00 (isothermal case) , we recover the dimension- 
less parameter 77 = kM/AirT of Sec. 14.7.11 On the other 
hand, in the homogeneous phase {B = 0), the energy iQ is 
given by E = ^ J pd9 = ^MT. Therefore, the normalized 
energy e = SirE/kM"^ is expressed in term of 77 according 
to e = i. In terms of these dimensionless parameters, the 
uniform phase is formally nonlinearly dynamically stable 
for 



V < Vcrit = 7i e > Cc 



1 

7 



(156) 



and linearly dynamically unstable otherwise. Note how the 
critical energy and temperature are simply expressed in 
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terms of the polytropic index 7. For n — s- +00, we recover 
the case of isothermal stellar systems with ecrit = 1 and 
Vcrit = 1- Note that the line {ecrit, Vcrit) coincides with the 
line i? = in Fig. |3 We thus clearly see how the series of 
equilibria for polytropic distributions places itself in the 
(CjTy) plane (we just have to displace the critical point 
(ccrit, Vcrit) along the line B = as sketched in Fig. I12|l. 

We also emphasize that, using the criterion H115|) we 
have obtained the condition of nonlinear dynamical insta- 
bility l|156(l for stellar polytropes with almost no calcula- 
tion. Of course, the same result can be obtained from the 
criterion (|116|l by explicitly performing the integral. The 
criterion of nonlinear dynamical stability H115(l that we 
have found is simpler, albeit equivalent. Moreover, it has 
a more physical interpretation since it is expressed as a 
condition on the velocity of sound in a gas with the same 
density distribution as the original kinetic system. 

For the polytropic distribution Ij 153(1 with uniform den- 
sity, the dielectric function can be written 



(157) 



where we have introduced the function 

a;[l 



/27r n 



in - 1/2) 



£l/2ln-3/2 

^ dx (158) 



with Wn{0) = 1. The range of integration is such that 
the term in brackets remains positive. For n — > +00, we 
recover the VF-function ((1^ . For w = 0, we obtain the 
critical temperature T = r„ as in a polytropic gas. As in 
Sec. 14.7.11 we look for solutions of the dispersion relation 
e(l, Lu) — in the form uj — iX where A is real. For T < r„, 
the system is unstable and the growth rate A > is given 

by 



T 



A 



where 
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Bn 
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fls/TT \ 2 



1 - 
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(159) 



dt 



For T > r„, the system is stable and the damping rate 
7 — —A > is given by 

T 



(160) 



where 



Fn{x) 



1 



G„(x)-i +i?„(x) 



with 



Rn{x) 
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a;2 1 










n + 1 



(161) 



(162) 



This additional term comes from the residue theorem when 
the pulsation cj = — 17 lies in the lower half of the complex 
plane. For n — > -l-oo, we recover the G and F functions 
((132|l and ((138|l . The dependence of the growth rate and 
decay rate with the temperature is shown in Figs. [T51 and 
1141 For n = 00, we recover the isothermal case of Fig. 1111 
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Fig. 13. Growth rate and decay rate for stellar polytropes 
and polytropic stars in the framework of the HMF model. The 
index goes from n = 1 to n = 100. We have also shown the case 
n = 1/2 (water-bag distribution). The critical temperature is 
smaller than for an isothermal gas (T„ < Tc). 
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Fig. 14. Growth rate and decay rate for stellar polytropes and 
polytropic stars in the framework of the HMF model. The in- 
dex goes from n — —1.5 to n = —100. The critical temperature 
is larger than for an isothermal gas {Tn > Tc). 



4.7.3 Fermi or water-bag distribution 

For n — 1/2, the distribution function 1(143(1 is a step 
function: /(e) = rjo if —vo{d) < v < vo{0) and /(e) = 
otherwise. This is similar to the Fermi distribution at 
T = describing cold white dwarf stars in astrophysics 
(Chandrasekhar 1939). This is also called the water-bag 
distribution in plasma physics (when vq is independent on 
0). The density and the pressure are given hy p = 2r]oVo 
and p = {2/3)r/oVQ. This leads to a polytropic equation 
of state p — Kp^ of index n = 1/2 and polytropic con- 
stant K = l/( 127^0 )• For a homogeneous system, we have 
the relation M = 47r77ofo- Then, combining the preceding 
relations, we find that the velocity of sound is Cs = wq- 
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Therefore, the system is formally nonlinearly dynamically 
stable if 



kM 
4^' 



(163) 



and linearly dynamically unstable otherwise. Noting that 
the kinetic temperature is T = Wq/S, we check that the 
above result returns (|98|l with 7 = 3. Thus, rjcrit — 3 and 
Ccrit = 1/3. Once again, these results have been obtained 
with almost no calculation. This is an advantage of for- 
mula (|115|l with respect to formula 1116|l . 

On the other hand, using f'{v) — riQ[5(v + wq) — 5{v — 
uq)], the dielectric function (|117|l can be written 

e(l,^) = 1 - ^W^i/2(^/VT), (164) 

with 



1 Q - 



(165) 



We look for solutions of the dispersion relation e(l, w) — 
in the form lo — Q where fl is real. This solution only 
exists for T > T1/2 and corresponds to an oscillatory so- 
lution Sf ~ e*^*. The pulsation is given by 



n = ±V3{T - = ±^vl - ^. (166) 

We now consider the case uj = iX where A is real. This 
solution only exists for T <Ti/2 and 



(167) 



The case A > corresponds to a growing (unstable) mode 
Sf ^ e'^* and the case A = —7 < corresponds to a 
damped mode Sf ~ e'"^*. We note that for this special 
case n = 1/2, the growth rate and the pulsation period of 
the stellar system are the same as for the corresponding 
barotropic gas, see Eq. l|S7|l . The results (|166|l and (|167|l 
have been previously derived by Choi & Choi (2003). They 
are recalled here for sake of completeness and because we 
will need them in Sec. 15.31 



5.1 The evolution of the whole system: the Landau 
equation 

There arc different methods to obtain a kinetic equation 
for the distribution function f(9,v,t). One possibility is 
to start from the N-body Liouville equation and use pro- 
jection operator technics. This method has been followed 
by Kandrup (1981) for stellar systems and by Chavanis 
(2001) for the point vortex gas. We shall first consider an 
application of this theory to the HMF model (Chavanis 
2003b). In the large N limit and neglecting collective ef- 
fects, the projection operator formalism leads to a kinetic 
equation of the form 

C f) Pi 

X .F(l->0,<-r)— +.F(0^ l,t-r) — 
ov ovi 

xf{ei,vut-T)f{e,v,t-T). (168) 

Here, f{9,v,t) = NPi{6,v,t) is the distribution func- 
tion, {F){9, t) is the (smooth) mean-field force and J-{\ 
0, t) = F(l ^ 0, t) - {F){e, t) is the fluctuating force cre- 
ated by particle 1 (located at 0i,fi) on particle (lo- 
cated at 0, v) at time t. Between t and t — t, the parti- 
cles are assumed to follow the trajectories determined by 
the slowly evolving mean-field {F){6,t). Equation H168|l 
is a non-Markovian integrodifferential equation. We insist 
on the fact that this equation is valid for an inhomoge- 
neous system while the kinetic equations presented below 
will only apply to homogeneous systems. Unfortunately, 
Eq. (jl68|l remains too complicated for practical purposes 
and we will have to make simplifications. If we consider 
a spatialy homogeneous system for which the distribution 
function / = f{v,t) depends only on the velocity, and if 
we implement a Markovian approximation, the foregoing 
equation reduces to 



dt dv 



r+00 

J dr J dOidviFil 0,t) 



5 Collisional relaxation of stellar systems 

The Vlasov equation H104|) can be obtained from the BBGKY, 
hierarchy, issued from the Liouville equation (0), by us- 
ing the mean-field approximation which is valid in the 
limit A'' +00 with 77 and e fixed. We would like now to 
take into account the effect of correlations between parti- 
cles in order to describe the "collisional" relaxation. We 
shall develop a kinetic theory which takes into account 
terms of order 1/A^ in the correlation function. 



xi^(l ^ 0, i - r) - f{v^,t)f{v, t), (169) 

where F(l ^ 0,t) = -^sm{0{t) - ei{t)). We thus need 
to calculate the memory function 

M = J dr J dOiFil 0, t)F{l -> 0, f - t) 

U2 r+00 r 

= ^J dr J d0iain{0-ei)sm[0{t-T)~9i{t-T)], 



(170) 

where 0i(t — r) is the position at time t — t oi the i-th 
particle located at 0i = 0i{t) at time t. Since the system is 
homogeneous, the mean force acting on a particle vanishes 
and the average equations of motion are 0{t — T) = 9 — vt 
and 0i{t ~ t) ^ 01 — viT. Thus, we get 

Af = — / dr d(p sin (j) sm{(t) - ut), (171) 
47'' Jo Jo 
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where 



9 — 6i and u 



vi . The integration yields Kubo formula 



t2 r + co .2 

M ^ — J dT cos(ur) = —S{u). (172) 
Therefore, the kinetic equation (|169|l becomes 

This equation can be considered as the counterpart of the 
Landau equation describing the "coUisional" evolution of 
stellar systems such as globular clusters (in that case, the 
system is not homogeneous but the collision term is often 
calculated by making a local approximation). The Lan- 
dau collision term can also be obtained from the BBGKY 
hierarchy at the order 0{1/N) by neglecting the cumu- 
lant of the three-body distribution function of order 
(see Chavanis 2004b). For the HMF model, and for one 
dimensional systems in general, we find that the Landau 
collision term vanishes. This is because the diffusion term 
(first term in the r.h.s.) is equally balanced by the friction 
term (second term in the r.h.s.). A similar cancellation 
of the collision term at order is found in the case of 
2D point vortices when the profile of angular velocity is 
monotonic (Dubin & O'Neil 1988, Chavanis 2001, Dubin 
2003). Therefore, after a phase of violent relaxation, the 
system can remain frozen in a stationary solution of the 
Vlasov equation for a very long time, larger than NtD. 
Only non-trivial three-body correlations can induce fur- 
ther evolution of the system. However, their effect is diffi- 
cult to estimate. Note that the collision term of order 1/iV 
may not cancel out in the case of inhomogeneous systems 
and for the multi-species HMF model (see Secdl. 



5.2 The evolution of a test particle in a thermal bath: 
the Fokker-Planck equation 

Equations H168|l and p69(l can also be used to describe 
the evolution of the distribution function P(u, t) of a test 
particle evolving in a bath of field particles with static 
distribution function /i(fi). In that case, we have to con- 
sider that the distribution function of the bath is given, 
i.e. /i = fo{vi). The evolution of P{v,t) is then governed 
by the equation 



dp_ _ d_ r+°° 

dt dv 



P + CO p 

J dr J d9idviF{l 0,t) 



xF{l -> 0,t-r) 



d d 

dv dvi 



foivi)P{v,t). (174) 



Equation Ij 174(1 can be written in the form of a Fokker- 
Planck equation 



dt dv \ dv 



(175) 



The two terms of this equation correspond to a diffusion 
and a friction. The diffusion coefficient is given by the 



D = 



+ 00 



dT{F{t)F{t - T))d7 



(176) 



The friction can be understood physically by developing 
a linear response theory. It arises from the response of 
the field particles to the perturbation induced by the test 
particle, as in a polarization process (see Kandrup 1983). 

Introducing the memory function H172|l . the Fokker- 
Planck equation (|174|l can be rewritten 

f = tI;/''"*--)(I;-^)a(-)^-(...'). 

(177) 

It can be put in the form H175|l where the coefficients 
of diffusion and friction are explicitly given by (Chavanis 
2003b): 

D=^ I dv,5{v - «i)/o(«i) = ^/o(«), (178) 



' 4 



dviS{v - vi)^{vi) 
dv 



-foiv). (179) 



More precisely, comparing Eq. ()175|l with the general Fokker- 
Planck equation 
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we find that 
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Using Eqs. H178|) and H179|) . we obtain 

1 (Av) 



dD 

dv 



''^-YAt 



(180) 



(181) 



(182) 



Therefore, rj represents half the friction force. This is the 
same result as in the case of Coulombian or Newtonian 
interactions (see, e.g., Chavanis 2004a). Using Eqs. H178|l 
and (|179|l . the Fokker-Planck equation H175() can be rewrit- 
ten 



dP 
'dt 



d_ 

dv 



p 



dv 



(183) 



with the initial condition P{v,t = 0) = S{v — Vq). It 
describes the evolution of a test particle in a potential 
U{v) = — ln/o(u) created by the field particles. A similar 
equation is found for 2D point vortices in position space 
where the friction is replaced by a drift (Chavanis 2001). 
When /o(w) is the Maxwellian H18|) . corresponding to a sta- 
tistical equilibrium state (thermal bath approximation), 
Eq. H183|l takes the form of the Kramers equation 
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as in the theory of Brownian motion (Risken 1989). We 
note in particular that the friction coefficient is given by 
the Einstein formula ^ — Df3. However, in the present 
context, the diffusion coefficient depends on the velocity 
and, in the ballistic approach that we have considered, is 
given by 



10 



1/2 



(185) 



We note that P{v,t) always converges to the distribution 
function of the bath NP{v,t) — > fo{v) for t — > +oo while 
for 31? self-gravitating system, this is the case only when 
/o is the statistical equilibrium distribution. 

Finally, again neglecting collective effects, a simple cal- 
culation shows that the temporal correlations of the force 
are 
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Fig. 15. Dependence of the diffusion coefficient D{v) on 
the velocity i; for different values of the inverse temperature 
7? = /3//3c. The normahzation constant is Do = ^(/3/27r)^''^ 
corresponding to D{0) for rj — 0. We note that D{0) increases 
as 77 increases and that it diverges at the critical point ri — 1. 



where D{t) is the Fourier transform of D{v). For the Maxwellian 



distribution, we get 



fo{'>^i)i we obtain Eq. H183|) where the expression of the 
diffusion coefficient is now given by 



{F(0)F(t)) = —cos(vt)e~^, 
An 



(187) 



D 



which seems to indicate a gaussian decay of the correla- 
tions. 



5.3 Collective effects: the Lenard-Balescu equation 

The kinetic theory developed previously (Chavanis 2003b), 
while useful as a first step, is however inaccurate because 
it is based on a ballistic approximation and ignores collec- 
tive effects which are non-negligible for the HMF model 
close to the critical temperature (Bouchet 2004). In the 
case of 3D stellar systems, collective effects have only a 
weak influence on the kinetic theory and they are often 
neglected. This implicitly assumes that the size of the sys- 
tem is smaller than the Jeans length (recall that the Jeans 
length plays the role of the critical temperature in the 
present context). In general, collective effects can be taken 
into account by developing a kinetic theory as in the case 
of plasmas (Ichimaru 1973). Noting that the HMF model 
is the one Fourier mode of a one dimensional plasma, Ina- 
gaki (1996) proposed to describe the colhsional evolution 
of the system by the corresponding form of the Lenard- 
Balescu equation. It can be written 



4 \e{l,vW 



(189) 



It differs from the preceding expression H178|l due to the 
occurence of the term |e(l,z;)p in the denominator which 
takes into account collective effects. For the Maxwellian 
distribution, the dielectric function is given by Eqs. H128|l 
and H130|) . This leads to the expression of the diffusion 
coefficient 



D{v) 
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l/2p-/3%- 



[1 - ijA{^vW + f ??2/5f2e-/3''' 



(190) 



with A{x) 1 - 
,.2 



dj_ _ }f_d_ r 6{v 

'dt~Td^j |e(i,w)P 



ov 



(188) 



where e(l, v) is the dielectric function (|117() . We note that 
the collision term again cancels out. However, if we use 
this equation to describe the evolution of a test particle in 
a thermal bath, as we did in Sec. 15.21 by replacing /i by 



^2 /p^e"2 du. We note that A{x) = 
1 — X"' + ... for a; — > and A{x) ~ —1/x'^ for x +00. 
Therefore, the diffusion coefficient behaves as Eq. (|185|) for 

V ^ +00 and tends to ^(/3/27r)i/V(l - tj)^ for v ^ 
and ry < 1. At the critical temperature 77 = 1 it diverges 

as D{v) ~ ^(/3/27r)i/V(/3w^) for v ^ 0. Its behaviour is 
represented in Fig. 1151 

The expression (I190|l of the diffusion coefficient was 
obtained by Bouchet (2004) in a different manner, by an- 
alyzing the stochastic process of equilibrium fluctuations. 
This approach was then generalized to an arbitrary dis- 
tribution function by Bouchet & Dauxois (2004), leading 
to equation (|189|l . In fact, formulae H189(l and l|190(l corre- 
spond to the one dimensional version of the diffusion coef- 
ficient computed in plasma physics (Ichimaru 1973) when 
the potential of interaction is truncated to one Fourier 
mode. The general expression of the diffusion coefficient 
and of the Fokker-Planck equation is given in Chavanis 
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Fig. 16. Same as Fig. 1151 for a polytropic distribution with 
index n = 1. Here rj = kM/AnT where T is the kinetic temper- 
ature. In this case, the critical point is rji = 2. The diffusion 
coefficient vanishes at the maximum velocity Vmax = V 4r. 
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Fig. 17. Same as Fig. 1151 for a polytropic distribution with 
index n = —2. Here r; = kM/AnT where T is the kinetic tem- 
perature. In this case, the critical point is rj-2 ~ 1/2. This 
figure is very similar to Fig. 1151 except that the diffusion coef- 
ficient decreases algebraically. 



(2004b) for an arbitrary potential of interaction in D di- 
mensions. 

We emphasize that the previous results are valid for 
an arbitrary distribution function fo{v) of the bath pro- 
vided that it is stable with respect to the Vlasov equation. 
This is because, as we shall see, the relaxation time of 
the "field particles" (bath) towards statistical equilibrium 
(Maxwellian) is longer than the relaxation time of a "test 
particle" towards /o, so that the distribution function /o 
of the bath can be considered as "frozen". The general 



expression of the diffusion coefficient can be written 



D{v) 



(191) 



where V stands for the principal value. Its asymptotic 
behaviour for v — s- +oo is always given by Eq. H178|l . As 
a complementary example to Fig. 1151 we have plotted in 
Figs. 1161 andllTlthe diffusion coefficient corresponding to a 
polytropic distribution of index n — \ and n = — 2 (in that 
case /o is given by Eq. H153|l where p and T are uniform) . 
Another example is provided by the water-bag model for 
which an explicit expression of D{v) can be given. Using 
Eqs. (|164|) . ()165|l and (|166|l . it can be written conveniently 
as 



(192) 



for V < vq and D — otherwise. The diffusion coeffi- 
cient diverges when the velocity of the particle v is in 
phase with the frequency < vq of the wave arising from 
the slightly perturbed distribution of field particles (see 



Sec. 14. 7. .3(1 . This divergence occurs because, for the wa- 
ter bag distribution, there exists purely oscillatory modes 
for a wide range of temperatures. In general, the diffu- 
sion coefficient ()189ll diverges only at the critical point for 
V — ujr with ,f'{ijJr) — 0; this precisely correspond to the 
criterion of marginal stability H120() . For example, for the 
gaussian distribution, we recover the divergence at T = Tf, 

for V = UJr = 0. 



5.4 The auto-correlation function 

We note that for high temperatures (i.e., 77 — > 0), Eq. 
(|19U|) reduces to the expression (|185|l found in the ballis- 
tic approach developed in Sec. 15. 21 This is because, in that 
limit, collective effects are weak with respect to the pure 
ballistic motion and e ~ 1 according to Eq. (|128|l . How- 
ever, the behaviour of the correlation function is different. 
Indeed, a direct analysis (Chavanis 2004b) shows that the 
temporal auto-correlation function of the force is given by 



(mm) 



47r 



cos[{v-vi)i\ 
12— /o(wi)awi 



= 2cos{vt)D{t). (193) 

Substituting this expression in the Kubo formula H176|l 
and using S{v — vi) ~ ^ J^^°° cos[(w — vi)t]dt, we recover 
Eq. (|189|l . In the ballistic approximation, the correlation 
function p87|l is gaussian while the exact treatment tak- 
ing into account collective effects shows that the decay 
of the fluctuations is in fact exponential with a decay 
exponent 7(/3) = {2/ /S^/'^F-'^iT]) where F{x) = 1/(1 + 
y/nxe^ erfc(— x)) is the i^- function (|138|l . This result was 
obtained by Bouchet (2004) by working out the integro- 
differential equation satisfied by the auto-correlation func- 
tion. We shall present here an alternative derivation (Cha- 
vanis 2004b) which will make a clear link with the decay 
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exponent appearing in the linear stability analysis of the 
Vlasov equation in Sec. 14.7.11 Noting that the correlation 
function is proportional to the Fourier transform of the 
diffusion coefficient, according to Eq. H19^^|l . we can ob- 
tain the expression of 7 by determining the pole of D(v) 
in Eq. H189() . Setting v — iX where A is real, we find after 
some calculations that 



|e|2(l,zA) = e(l,zA)e(l,-zA) 



(194) 



where we recall that e(l, iX) is given by Eq. Hl^Ufl . Clearly, 
|ep(l,zA) is an even function of A. We need to determine 
the values of A for which this function vanishes. Since 
77 < 1, we find that A = ±7 where 7 > is determined by 
e(l, —17) = 0. Therefore, 7 is the damping rate of the sta- 
ble perturbed solutions of the Vlasov equation; it is given 
by Eq. (|137|l . Next, we consider A = ±7 -|- e. Expanding 
Eq. H194|) for e ^ 1, we find after elementary calculations 
that 



D{v) 



y2 _|_ ^2 ■ 



2T 



1 



(195) 



for V — > ±17. Therefore, for t — > +00, the correlation 
function H193() is the Fourier transform of a Lorentzian so 
it decays like 



{F{0)F{t)) 



1 



7 |F'( 



^)l 



%/2T 



with 



7=(2//3)i/2F-i(,y). 



cos{vt)e-^\ (196) 



(197) 



The exponential decay of the correlation function corre- 
sponds to the Markovian limit of the stochastic process, 
thereby justifying the Markovian approximation in the ki- 
netic theory. This is quite different from the correlations of 
the gravitational force which decay as t^^ (Chandrasekhar 
1944). This slow decay may throw doubts on the validity 
of the ordinary Landau equation, based on a Markovian 
approximation, used to describe stellar systems (see Kan- 
drup 1981 for a detailed discussion). For the HMF model, 
the decay exponent 7(/3) depends on the temperature. It 
diverges like 7 ~ \/2T\nT as T ^ +00. This is why the 
correlation function is gaussian in the treatement neglect- 
ing collective effects (see Sec. I5.2|l . On the other hand, 
7 ~ (8/fcM)i/2(T-rc) as T -> T+, so that the correlation 
function decreases very slowly close to the critical temper- 
ature. This may invalidate the Markovian approximation 
close to the critical point. The Fokker-Planck equation 
H184(l with the diffusion coefficient l|19()(l has been recently 
investigated by Bouchet & Dauxois (2004). In this paper, 
using the rapid decay for large v of the diffusion coeffi- 
cient H190(l . the numerically observed (Yamaguchi 2003, 
Pluchino et al 2004) anomalous algebraic decay of the ve- 
locity autocorrelation function is explained and algebraic 
exponents are explicitely computed. For a large class of 
bath distribution function /o, this may also lead to anoma- 
lous diffusion of angles 9. 



5.5 The relaxation timescale 

Let us consider the relaxation of a test particle in a ther- 
mal bath. Due to the rapid decrease of D{v) for large v, the 
spectrum of the Fokker-Planck equation (|184|l has no gap, 
and it exists no exponential relaxation time (see Bouchet 
& Dauxois 2004 for a detailed discussion). A time scale 
will however describe relative relaxation speeds for differ- 
ent values of the temperature T. We shall obtain an esti- 
mate of this time scale. Ignoring collective effects in a first 
step, this process is described by Eqs. (|184(l - (|185|l . If the 
diffusion coefficient were constant, we would deduce that 
the dispersion of the particles increases as {{Av)"^) = 2Dt. 
Introducing the r.m.s velocity Vm — (w^)^/^, we define the 
relaxation timescale tr such that ((Z\u)^) = vf^. This leads 
to tr — v'^/2D. Since D depends on v, the description of 
the diffusion process is more complex. However, the for- 
mula 



2D{v„^)' 



(198) 



provides a useful estimate of the speed of relaxation. For 
the Maxwellian distribution for which = 1 / we get 

(199) 



0.121pfc2 



We can also estimate the relaxation timescale by = 
1/^, where ^ is the friction coefficient. Using the Einstein 
relation ^ = Df], this yields = 2t,.. Finally, setting w — 
v/{^/2vm), we can rewrite Eq. H184() in the dimensionless 
form 



dP _ J__9_ 

dt tji dw 

where G{w) ~ e' 



and 



tR 



0.05pfc2' 



(200) 



(201) 



which provides a useful time renormalization factor in the 
Fokker-Planck equation. If we use expression (|190|l of the 
diffusion coefficient, we need to multiply the relaxation 
timescale by |e(l,'i;,„)p ^ {1 + 2.2bj^f -f 0.57877^. In par- 
ticular, close to the critical temperature (i.e. = 1), the 
relaxation timescale is multiplied by ~ 11.1. 

If we now introduce a dynamical timescale through tlie 
relation 



27r 

to — — , 



(202) 



we find that tn/tD ^ v^/Nk^ - l/Nk'^0^ ~ N/rf ~ N 
in the thermodynamic limit of Sec. 12. II Therefore, the re- 
laxation time of a test particle in a thermal bath increases 
linearly with the number N of field particles. More pre- 
cisely, we have 



tR 



0.127 

9 



Nt 



D- 



(203) 
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The collision term in Eq. 12()()|) is of order in the 
thermodynamic limit N +00 with 77 fixed. This scaling 
also holds for the collision term in the Landau equation 
1)173(1 noting that k 1/N and / ~ if we take ~ 1- 
It represents therefore the first correction to the Vlasov 
limit in an expansion in 1/iV of the correlation function. 
However, contrary to the relaxation of a single particle in 
a thermal bath, the relaxation time of the whole system 
is not treiax ^ Nt D bccausc, as we have seen, the collision 
term cancels out at the order Therefore, the relax- 

ation time of the whole system is larger than Ntjj. This 
is consistent with the finding of Yamaguchi et al. (2004) 
who numerically obtain treiax ^ N^''^to- 

It is interesting to compare this result with other sys- 
tems with long-range interactions. For stellar systems, the 



Chandrasekhar relaxation time scales as t^ 



N 



to. 



l-relax ■ - in JV ' 

It corresponds to the Nt d scaling polluted by logarithmic 
corrections. This is the relaxation time of a test star im- 
mersed in a bath of field stars as well as the relaxation time 
of the whole cluster itself (in the absence of gravothermal 
catastrophe). For the point vortex gas, the collision term 
in the kinetic equation cancels out at the order \/N when 
the profile of angular velocity is monotonic. Therefore, the 
relaxation time of the whole system is larger than Nt^. In 
fact, it is not clear whether the point vortex gas ever re- 
laxes towards statistical equilibrium (Chavanis 2002a, Du- 
bin 2003). By contrast, the relaxation time of a test vortex 



In AT 



to 



in a thermal bath of field vortices scales as treiax 
(Chavanis 2002a). 

6 Self-attracting Brownian particles: the 
BMF model 



The Hamilton equations (^3) describe an isolated system 
of particles with long-range interactions. Since energy is 
conserved, the fundamental statistical description of the 
system is based on the microcanonical ensemble. It can be 
of interest to consider in parallel the case where the system 
is stochastically forced by an external medium. We thus 
introduce a system of Brownian particles with long-range 
interactions which is the canonical version of the Hamilto- 
nian system This could be called the BMF (Brownian 
Mean Field) model. Likewise in the case of 3D Newto- 
nian interactions, a system of self-gravitating Brownian 
particles has been recently introduced and studied (see 
Chavanis et al. 2002 and subsequent papers). 



6.1 Non-local Kramers and Smoluchowski equations 

We consider a one-dimensional system of self-attracting 
Brownian particles with cosine interaction whose dynam- 
ics is governed by the A^-coupled stochastic equations 




where — is a friction force and Ri{t) is a white noise 
satisfying 

{R,{t)) = 0, {Rr{t)Rj{t')) = 5,,5{t - t'), (205) 

where i — 1, ...,N refer to the particles. The particles in- 
teract through the potential U{9i, ...,6n) — J2i<j^i^i ~ 
9j) where u{9i — 6j) = — ^ cos(6'i — 9j). We define the in- 
verse temperature (3 = 1/T through the Einstein relation 
^ = Dp. The stochastic model (I204ll - (|205ll is analogous 
to the model of self-gravitating Brownian particles intro- 
duced by Chavanis et al. (2002). For this system, the rele- 
vant ensemble is the canonical ensemble where the temper- 
ature measures the strength of the stochastic force. The 
evolution of the TV-body distribution function is governed 
by the A^-body Fokker-Planck equation 



N 



N 



dt 



El OPn , ^ OPn 
I Vi + Fi 



1=1 

N 



39, 



dv,. 



dvi 



dv. 



(206) 



where Fi — — Jf-- The stationary solution corresponds to 
the canonical distribution 



Pn 



(207) 



We note that the canonical distribution (|207|l is the only 
stationary solution of the A^-body Fokker-Planck equation 
while the microcanonical distribution (jS)) is just a partic- 
ular stationary solution of the Liouville equation (see Sec. 

For the system pHjl - lpUS)) . the equilibr ium canonical 
distribution does not rely, therefore, on a postulate. In the 
thermodynamic limit N — s- -l-oo with 77 = j3kM /47r fixed, 
one can prove that the N particules distribution function 
factorizes and that the mean-field approximation is exact 
(Chavanis et al. 2004b, Chavanis 2004b). The evolution of 
the one particle distribution function is then governed by 
the non-local Kramers equation 



'dt^'"'d9~'d9d^~lh\ d^^^-'^J' 



(208) 



which has to be solved in conjunction with Eq. (jSJ. These 
equations have been considered, independently, by Choi 
& Choi (2003). 

To simplify the problem further, we shall consider a 
strong friction limit ^ — > +00 or, equivalently, a long time 
limit t ^ In that case, we can neglect the inertia of 
the particles and the stochastic equations H204|l reduce to 

d9 

' '^^—Ui9l,...,9N) + V^*R^{t), (209) 



dt 



39, 



where ^ — 1/^ is the mobility and £>* = 7?/^^ = T/^ is 
the diffusion coefficient in physical space. The evolution 
of the iV-body distribution function is governed by the 
A^-body Fokker-Planck equation 



dPj 



N 



2DR,{t), (204) -qT 



y— 



(210) 
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The stationary solution corresponds to the canonical dis- 
tribution in configuration space 



P, 



1 



N 



-/3c/(ei,...,ejv) 



(211) 



In the thermodynamic limit N — )■ +oo with rj ~ [3kM /At: 
fixed, the mean-field approximation is exact and the evo- 
lution of the one particle distribution function is governed 
by the non-local Smoluchowski equation 



dp d 



(212) 



where <P is given by Eq. ©. Alternatively, the Smolu- 
chowski equation (|212|l can be obtained from the non-local 
Kramers equation (|208|l by using a Chapman-Enskog ex- 
pansion in power of 1/^ (Chavanis ct al. 2004a). In that 
approximation, the distribution function is close to the 
Maxwellian 



fiO,v,t) 



1 



V2^ 



p{d,t)e- 



(213) 



and the evolution of the density is governed by the non- 
local Smoluchowski equation. The equations H2U8|) and 
H212(l conserve mass and decrease the Boltzmann free en- 
ergy, i.e. Fb < 0. This is the canonical version of the 
H-theorem (Chavanis 2003c,2004a). 



Chavanis et al. (2002) for self-gravitating Brownian parti- 
cles. We note that the eigenvalue equations (|91(l and H215|l 
for the Euler model and the Brownian model only differ 
by the substitution X^. Therefore, the results of Sec. 

I3.5.1l can be directly extended to the present context. 

For the uniform phase, the most destabilizing mode 
(n = 1) is 



Sp = ai cos 9 e 



At 



with a growth rate 



(216) 



(217) 



For T < Tc, the perturbation grows expontentially while 
it is damped exponentially for T > Tc. 

Considering now the clustered phase and using a per- 
turbative approach similar to that of Appendix IXIfor T — > 

(not detailed), we find that the perturbation is damped 
exponentially with a rate 



(218) 



From Eqs. H217(l and 1)218(1 . we find that, starting from 
a homogeneous solution at T = T~ , there is first an ex- 
ponential growth on a timescale ^{Tc — T)^^ . Then, non- 
linear terms come into play and the system relaxes towards 
a clustered state on a timescale (l/2)^(rc — T)~^ . 



6.2 Linear stability 

The stationary solutions of Eq. (|212|l are given by Eq. 
((19(1 . They extremize the Boltzmann free energy Fb = E— 
TSb with and at fixed mass and temperature. 
Furthermore, only stationary solutions that minimize the 
free energy are linearly dynamically stable with respect to 
the non-local Smoluchowski equation (Chavanis 2003c). 
Therefore, thermodynamical and dynamical stability are 
clearly connected: the stable stationary solutions of Eq. 
((212(1 correspond to the canonical statistical equilibrium 
states in the mean-field approximation. 

Considering a small perturbation 6p around a station- 
ary solution of Eq. ((212(1 . we get 



dSp d 



T 



d6p 
80 



d<P dS<I> 



Writing Sp ~ e'*'* and introducing the notation l(50() . we 
obtain the eigenvalue equation 

d fldq\ k '^"^ 



de\pdej 2ttT 



q{e') cos{e - 9')de' = ^q, 

Tp 



(215) 



which is similar to the eigenvalue equations obtained in the 
preceding sections. These eigenvalue equations all coincide 
at the point of marginal stability A = implying that the 
onset of the instability is the same in all the models con- 
sidered. Equation ((215(1 is the counterpart of Eq. (27) of 



6.3 Local Fokker-Planck equation 

Before studying the non-local Smoluchowski equation ((212(1 - 
(|Sl numerically, it may be of interest to study a simpli- 
fied problem where is replaced by its equilibrium value 
= BcosO, where B is fixed. Considering again the 
strong friction limit, we obtain the classical (local) Smolu- 
chowski equation 



dp _ d 
dt ^09 



(219) 



This Fokker-Planck equation also appears in the model of 
rotation of dipoles in a constant electric field developed by 
Debye. The stationary solutions of this equation are given 
(214) by Eq. 1(19(1 . Considering a perturbation Sp ~ e"^* around 
a stationary solution and setting Sp — dq/d9, we obtain 
the eigenvalue equation 



d_ 

d9 



1 dq 

'pde 



AC 
Tp 



(220) 



This equation has the form of a Sturm-Liouville problem. 
For _B ^ 1, so that p is approximately uniform, we find 
to leading order that the eigenvalues are A„ = —n^T/^ 
where n = 1,2,... Using a procedure similar to that of 
Appendix ^ (not detailed), the next order correction is 
-Xn^/T = n^ + n^/[2{An^-l)]x'^ + ... where a: = l3B < 1. 
The eigenvalues of the Sturm-Liouville equation 1(220(1 are 
evaluated numerically in Fig. 1181 
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Fig. 18. Eigenvalues of the Sturm- Liouville equation 122UII . 



For B <^ 1, the Fokker-Planck equation (|219|l can be 
solved analytically. Assuming that, initially, the density is 
uniform, we find 



M 
2^ 



1 



1)(3B cose 



(221) 



For t ^ treiax — (.P (relaxation time), the density reaches 
its equilibrium value (|19|l . For short times i ^ 1, we have 



M f 1 

— 1 - -tB cose 

27rV e 



(222) 



Starting from a homogeneous solution, there is first a lin- 
ear growth followed by an exponential relaxation towards 
the clustered state on a timescale ^ i/T. This is illus- 
trated numerically in Fig. 1261 and compared with the time 
evolution of the non-local Fokker-Planck equation H212(l . 
The evolution is of course quite different since, as we have 
seen, the non-local Fokker-Planck equation displays an ex- 
ponential growth on a timescale ^(Tc — T)^^ followed by 
an exponential relaxation on a timescale (l/2)^(rc — T)^^. 
These timescales diverge as we approach the critical tem- 
perature Tc while there is no critical temperature for the 
local Fokker-Planck equation where B is fixed. The nu- 
merical solution of Eq. (|219|l is shown in Fig ^] and the 
analytical solution (|221|l in Fig.lJOlfor two different values 
of rj. 

We also note that the change of variables 



(223) 




Fig. 19. Evolution of the density profile according to the local 
Fokker-Planck equation for -q — 2 (corresponding to a; = 3.33) 
starting from a homogeneous solution (numerical simulation). 




Fig. 20. Evolution of the density profile according to the local 
Fokker-Planck equation for rj = 1.005 (corresponding to a; = 
0.2) starting from a homogeneous solution. In this case, we have 
used the analytical solution 122H which is valid for B <^ 1. 



In our case 



v{e) = -- 



Bcose+^B"^ sin^i 



(226) 



transforms the Fokker-Planck equation H219|) into a Schrodinger- 



like equation (with imaginary time) 



with 



ViO) ^ 2 



(224) 



(225) 



We can use this formalism to study the low temperature 
regime T — > 0. In that case, the equilibrium state is close 
to a Dirac peak centered on = tt so that we can expand 
the potential to leading order in a; = 6* — tt. Eq. H219|l 
becomes a Kramers-like equation 



dp d 
dt dx 



(227) 
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Fig. 21. Characteristics of Eq. For t +00, all the 

particles converge at ^ = tt and a Dirac peak is formed. 



Fig. 23. Characteristics of the damped Euler equation 123311 
for p = and <P = BcosO. The ratio ^/B = 1/2. 



20 - 
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Fig. 22. Evolution of the density profile at T = starting 
from a homogeneous distribution. 



The corresponding Schrodinger-like equation H224() is that 
of a harmonic oscillator 



dt 



T 



1 



1 



(228) 



Therefore, for T —t Q the eigenvalues are given by A„ = 
—nBjt,, i.e. \nS,/T = —nx (with n = 1,2, ...) and the cor- 
responding eigenfunctions are the Hermite polynomials. 

It is also possible to solve the problem exactly at T = 
by using the method of characteristics. Indeed, the Fokker- 
Planck equation becomes 



(229) 



istic and the equation of motion of a particle is 
dt £_ 

This equation of motion is readily solved and we get 



(230) 



9{t) = 2Arctan 



tan I — 



(231) 



where 9^ is the initial position of the particle. The car- 
acteristics are shown in Fig. [2] For t +00, all the 
particles converge at 6 = n and a Dirac peak is formed. 
The density profile is determined by the condition pqcIOq = 
p{e,t)de yielding p{9,t) = po/{d9/d9o). From Eq. P?T1) . 
we get 



Pi9,t) 



Po[l + tan^(6'/2)]e-»*/^ 
l + tan2(0/2)e-2St/? ' 



(232) 



This profile is shown in Fig. [53 at different times. 

To make the link between Figs. |^ and 1101 we can 

consider an intermediate hydrodynamical equation 



du 
'dt 



du 
'd9 



1 dp 
])d9 



'89 



(233) 



This could be called the damped Euler equation (see Cha- 
vanis 2003c). For ^ = we recover the Euler equation 
(|67|l and for ^ — )■ +00, using the continuity equation IjGGI) . 
we recover the Smoluchowski equation H212(l . For p = 
and <P = B cos 9, Eq. H233fl can be solved by the method 
of characteristics and the results are reported in Fig. 
This clearly shows the passage from Fig. ^3 to Fig. 1211 as 
the friction parameter ^ increases. 



6.4 Dynamics of Brownian particles in interaction 



which is equivalent to an advcction equation by a velocity 
field V ~ (B/^) sin 9. Therefore, the evolution is determin- 



We now turn to the evolution of a system of Brownian par- 
ticles in interaction described by the iV-coupled stochastic 
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Fig. 24. Evolution of the density profile according to the non- Fig. 26. Evolution of the different modes an(t) for the non- 
local Fokker-Planck equation. For 77 = 0.769 < 1, the homoge- local (full lines) and local (dashed lines) Fokker-Planck equa- 



neous solution is stable. 



tions. The control parameters are i) — 2^ x — 3.33. 




Fig. 25. For 77 = 2 > 1, the homogeneous solution (disk- 
like) is unstable and the system forms a cluster (bar-like) . The 
evolution is longer than with the local Fokker-Planck equation 
represented in Fig. 1191 as also shown in the next figure and 
explained in the text. 



equations H2U4|I . Despite all the simplifications introduced, 
this model remains a non-trivial and interesting model 
exhibiting a process of self-organization. In particular, it 
shows the passage from a homogeneous phase (disk-like) 
to a clustered phase (bar-like) under the influence of long- 
range interactions. Interestingly, the large N limit of this 
system is exactly described by an explicit kinetic equation 
(|2U8|) reducing to Eq. (|212|l for large times. By contrast, 
the kinetic equation describing the evolution of the HMF 
model towards statistical equilibrium is not known as 
the collision term cancels out at the order 1/iV within the 
approximations usually considered (see SecEjl- 



For the Brownian gas, we have to solve the integro- 
differential equation 



i^ = T^ + 7r^\pl sm{e -0')p{9\t)d9' 



27r 



Substituting the decomposition 

-\-co 

p = ao(t) + ^ a„(i) cos{n0) 



in Eq. lf2M|l . we find that 



ao 



M 
2^' 



^dai ^ k f M 02 

E V 1 a^ = — oi 

^ dt 2 H 27r 2 



. dan 
'dF 



-nai(a„_i - a„+i), 



(234) 
(235) 

(236) 
(237) 
(n> 2). (238) 



We note that the coefficient ai{t) is related to the magne- 
tization B{t), defined in Eq. (|10|l . by the formula B{t) = 
-|ai(t). According to Eqs. |(2U and 1221), the coef- 

ficients a„ are given at equilibrium by 



(239) 



Therefore, the static equations H237|) and H238|) with d/dt = 
coincide with the recursive relations satisfied by the 
Bessel functions /„(x). 

For T TjT, the coefficients scale B" with 

B 1. Therefore, after a transient regime of order l/Tn^, 
Eq. H238|l can be simplified in 



k 1 
4Tn 



aia„_i. 



(240) 
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In particular, for n = 2, we get 02 = ^a^. This shows 
that the second mode is slaved to the first. Substituting 
this in Eq. ((^Tfl) . we get 



dai fc2 3 

c.—, h (J — J-c)ai = ^fti • 

^ dt ^ ' 32T ^ 



(241) 



This equation is readily solved (it may be convenient to 
use p = as a variable) and we obtain 



where 



alit) 



A 



^(T, - T)e 



2(Te-T)t/5 



1 _|_ MLp2(T^-T)t/i 
1 -I- ^ 



fc^ai(0)3 
32T 



(242) 



(243) 



For t 



l-oo, we get ai(oo)^ = '^^{Tc — T), or equiva- 
lently = 8T{Tc — T). We thus recover the equilibrium 
result H32|l valid close to the critical point. The approach 
to equilibrium for t — )■ -\-oo is governed by 5ai (t) / ai {00) ~ 
^i67^g-2(T,-T)t/? yielding the damping rate ^T^. On the 
other hand, for t <C treiax, one has ai{t) = ai(0)e(-^'^~-^)*/^ 
yielding the growth rate H217|l . Finally, at T = Tc, Eq. 
|P1T|| leads to 



ai(<) = ± 



1 



ai(0)2 



(244) 



so that the magnetization B{t) 



}ai{t) tends to zero 



algebraically as t for t +00. Using Eq. (|240|l . the 
coefiicients a„(i) are expressed in terms of ai(t) by 



k 

4T 



1 



(245) 



In the previous calculations, we have assumed for simplic- 
ity that the density profile is symmetrical with respect 
to the a;-axis. The general case is treated in Appendix IfI 
Away from the critical point, the non-local Smoluchowski 
equation has to be solved numerically. Some numerical 
simulations are shown in Figs. I24I26I Note also that for 
T ^ and for sufficiently large times, the density is 
peaked around 9 — tt and Eq. H234|) becomes equivalent 
to the local Kramers equation (|227() with B = kM /2-k. 
Therefore, for T ^ 0, the eigenvalues of Eq. (|215|l tend 
to A„ = —nkM/2TT^ (harmonic oscillator). Substituting 
^ A^, we deduce that, for T — > 0, the eigenvalues of 
Eq. (|91|l tend to Xf^ = —nkM/2'K (they are represented in 

Fig. El). 



7 The multi-species HMF model 

We finally briefiy comment on the generalization of the 
preceding results to the case of the multi-species HMF 
model. We thus return to the Hamiltonian equations ^ 
and account for the possibility of having particles with 



different masses. Stellar systems also possess a mass spec- 
trum, so this generalization has a counterpart in astro- 
physics. 

Considering first the statistical equilibrium state, a 
straightforward generalization of the counting analysis of 
Sec. O yields 



_ _ .."in 



(246) 



for the probability of the state {n^a}, where riia gives the 
number of particles with mass nia in the i-th macrocell. 
Therefore, the entropy S = \nW oi the multi-species gas 



^ J ma ma 



(247) 



where fa{0,v)d9dv gives the total mass of particles of 
species a in {0,v). The distribution function of the whole 
assembly is 



f{e,v) = Y,ja{e,v). 



(248) 



The statistical equilibrium state is obtained by maximiz- 
ing the entropy (|247|1 while conserving the total energy E 
and the mass Ma of each species of particles. This yields 



(249) 



which generalizes Eq. H18|l. Note that the inverse temper- 
ature (3 is the same for all species in accordance with the 
theorem of equipartition of energy. This clearly leads to 
a mass seggregation since the r.m.s. velocity of species a 



decreases with the mass: 
Eq. lf^I!)|l imphes 



T/ma- More precisely. 



/a(e) = Cafc[/fc(e)]"»/"% 



(250) 



where Cab is a constant. Therefore, heavy particles will 
have the tendency to occupy regions of low energy. Re- 
call, by contrast, that there is no mass seggregation in 
Lynden-Bell's statistical theory of violent relaxation for 
coUisionless systems (see Sec. 0J since the mass of the 
particles does not appear in the Vlasov equation (Lynden- 
Bell 1967). It would be of interest to study these problems 
of mass seggregation with the HMF model for which nu- 
merical simulations are simpler than with gravitational 
systems. 

Considering now the coUisional relaxation, a straight- 
forward generalization of the Lenard-Balescu equation to 
a multi-species system yields 



dfa ^ ^d_Y^ 



with 



4 dv 



dv 



,S{v-v') 
WW 



„dfa 

dv 



e(l,w) = l + 



k X ^ 

2 ^ 

b 



dv. 



(251) 



(252) 
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We now see that the coUision term does not vanish any- 
more when there are at least two different species. The 
diffusion and the friction experienced by a particle of one 
species are caused by collisions with particles of another 
species. Equation H251|) can be rewritten in the suggestive 
form 



dfa 

dt 



-y 

dv ^ 

b 



J-^ab-^ Uab-^ 

OV OV 



where 



Dab = 



mbfb 



4 Hl,vW 



(253) 



(254) 



where Kab is a constant independent on v. This equation 
is similar to Eq. H25U|I but, here, fa{v) is not necessarily 
the Maxwellian. Indeed, Eq. H259I) is satisfied by any dis- 
tribution of the form /a(w) = Aaexp[— /3max(?;)], where 
xiv) is determined by the initial conditions. 

We can use these results to study the relaxation of a 
test particle of mass m in a bath of field particles with 
mass nif. Neglecting collective effects for simplicity, we 
find that the equivalent of the Fokker-Planck equation 
(| 184(1 is now 



is the diffusion coefficient for species a due to collisions 
with species b and 



Dab = 



ma fa 



4 \e{l,v)r 



(255) 



is an "off-diagonal" diffusion coefficient. It corresponds to 
a friction force 



Va 



Dab dfb 



fa dv 

These diffusion coefficients satisfy the relation 

Dab 



Dab ^ 

ma fa mbfb ' 



(256) 



(257) 



Ja 



E 



' dfa dfb 

^ab-^ ^ab-^ 

OV OV 



= 0, 



(258) 



implying the following equilibrium relation between the 
densities 



dP 
'dt 



d_ 

dv 



f3mPv 



(260) 



with 



D{v) = r 



(3m j 



1/2 



(261) 



where n is the number density of field particles. The equi- 
librium distribution of the test particle is 



The timescale of collisional relaxation is 

,,3 



These results are similar to those obtained by Dubin (2003) 
for the multi-components point vortex gas in two dimen- 
sions. When the profile of angular velocity is non-monotonic, 
the vorticity profile evolves under the effect of long-range 
collisions caused by a process of resonance (Dubin & O'Neil 
1988, Chavanis 2002). When the profile of angular velocity 
is monotonic, there is no evolution for the single-species 
system. An evolution is, however, possible for the distri- 
bution function of each species in the multi-components 
system. In a sense, at the order l/N, the kinetic theory of 
the point vortex gas (evolution of the single-species sys- 
tem when the angular velocity profile is non-monotonic) is 
intermediate between the kinetic theory of stellar systems 
(evolution of the single-species system in any case) and 
the kinetic theory of the HMF model (no evolution of the 
single-species system). 

Therefore, the results of Dubin (2003) can be directly 
transposed to the present context. In particular, we note 
that the total distribution function X]a/a(^'^) ~ 
stationary so that the conservation of energy is trivially 
satisfied. On the other hand, a 7J-theorem can be proved 
for the entropy (PTzjl . i.e. S > 0. The equality S = 
corresponds to vanishing currents 



mf 



0.121nm2fc2' 



(262) 



(263) 



and tj, = 2{mf /m)tr. In dimensionless form, Eq. H260() can 
be rewritten 



dP _ 1 d 
dt tji dw 



G{w)[^ + 2—Pw 
dw mf 



with 



^ 

O.Obnm'jk'^ 



(264) 



(265) 



and G{w) = e . If we properly account for collective 
effects, the diffusion coefficient is given by 



Div) = 



where rj 



[1 - rjAiy^^^v)]'^ + f ?72^m/u2e-'3"/"' 



(266) 



fcM m f 
iirT 



On the other hand, for a distribution of 



the bath of the form fo{v), the Fokker-Planck equation is 

(267) 



dP _ d 
dt dv 



Div)[^~^P^ 
dv mf dv 



faiv) = Kab[fb{v)r-^"'\ 



with D{v) = ■^m//o(w)/|e(l, w)p. We note that the equi- 



(259) librium distribution of the test particle is Peq{v) = A/o(w)'"/™^ . 
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8 Conclusion 

In this paper, we have given an exhaustive description 
of the HMF model that recently appeared in statistical 
mechanics as a simple model with long-range interactions 
similar to self- gravitating systems. The originality of our 
approach is to offer an overview of the subject and to see 
how different models (Hamiltonian, Brownian, fluids,...) 
are related to each other. Other studies concentrate in 
general on a specific aspect of the problem. We think that 
putting all the models in parallel is illuminating because 
they are closely connected to each other so that a unified 
(and aesthetic) description can be given. These connec- 
tions were previously noted by one of us (P.H.C) in the 
case of 3D self-gravitating systems and it was natural to 
extend these results to the HMF model. A more general 
approach is given in Chavanis (2004b) for an arbitrary po- 
tential of interaction in D dimensions. The present paper 
can be seen as a particular application of this general for- 
malism for a one dimensional potential truncated to one 
mode. The main interest of the HMF model in this context 
is to yield simple explicit results. 

Another originality of our approach is to emphasize the 
connection between the HMF model and self-gravitating 
systems (and 2D vortices) although this link is only sketched 
in other papers, except in the early work of Inagaki. Many 
concepts and technics that are well-known in astrophysics 
have been rediscovered for the HMF model, sometimes 
with a different point of view. This is true in particu- 
lar for the notion of violent relaxation and metaequilib- 
rium states. In statistical mechanics, this has been ap- 
proached via a notion of "generalized thermodynamics" 
(Tsallis 1988) although it was understood early in as- 
trophysics (Lynden-Bell 1967, Tremaine et al. 1986) that 
these metaequilibrium states correspond to particular sta- 
tionary solutions of the Vlasov equations on a coarse- 
grained scale (Chavanis et al. 1996, Chavanis 2003a). Thus, 
our dynamical interpretation of Tsallis functional as a par- 
ticular H-function differs from the thermodynamical inter- 
pretation given by Boghosian (1996), Latora et al. (2002) 
and Taruya & Sakagami (2003). 

Concerning the interest of the HMF model for astro- 
physicists, we have shown that it exhibits the same types 
of behaviors as 3D self-gravitating systems while being 
much simpler to study because it is one dimensional and 
avoids complicated problems posed by the divergence of 
the gravitational potential at short distances and the ab- 
sence of a large-scale confinement. Thus, it distinguishes 
what is common to long-range interactions and what is 
specific to gravity. This comparative study should bring 
new light in the statistical mechanics of self-gravitating 
systems which has long been a controversial subject. Other 
simplified models of gravity have been introduced such as 
the parallel planar sheets of Camm (1950), the concen- 
tric spherical shells of Henon (1971) or the toy models 
of Lynden-Bell & Lynden-Bell (1977) and Padmanabhan 
(1990). These toy models have often allowed advances in 
the description of more realistic self-gravitating systems 
that are difficult to study in full detail. We think that. 



similarly, the HMF model should find its place in the as- 
trophysical literature. 

Another interest of the HMF model is to allow to study 
in great detail what happens close to the critical point. In 
the HMF model, the potential is truncated to one mode 
n = 1 and there exists a critical temperature Tc = kM/Air 
below which the uniform phase is unstable. For infinite 
homogeneous self-gravitating systems, there is a continu- 
ous spectrum of modes but there exists a critical wave- 
length A J — {nT/GmpY^^ (depending on the tempera- 
ture) above which the system is unstable. Alternatively, 
if we fix the size of the system (for example by placing it 
in a box of radius R) , the maximum wavelength is R and 
the Jeans instability criterion now determines a critical 
temperature Tc GmM / {ttR) below which the system 
becomes unstable (we have used Xj — R and p ~ M/R^). 
In fact, in the case of box-confined gravitational systems, 
we must consider that the gaseous phase is inhomogeneous 
and use a more precise stability analysis (Chavanis 2002b) 
yielding Tc = GmM / {2.b2R). These remarks show that 
the critical temperature in the HMF model plays exactly 
the same role as the critical temperature in finite isother- 
mal spheres. Now, in the framework of the HMF model, it 
is possible to study how the mean-field results are altered 
close to the critical point due to collective effects. This 
is more complicated for self-gravitating systems because 
they are inhomogeneous. However, on a qualitative point 
of view, we expect similar results: divergence of the two- 
point correlation function like in Eq. (|61|l . divergence of 
the force auto-correlation function like in Eq. (|64|l . alter- 
ation of the diffusion coefficient and increase of the relax- 
ation time like in Eq. 11901) , increase of the decorrelation 
time like in Eq. 1197|) ... These problems are difficult to 
study for self-gravitating systems but they are of consider- 
able importance. They have never been discussed in detail 
because it is usually implicitly assumed that the system 
size is much smaller than the Jeans scale (or the temper- 
ature much larger than Tc) so that collective effects are 
neglected and the dielectric function is approximated as 
e ~ 1. The present simplified study is a first step to un- 
derstand the failure of the mean-field approximation close 
to the critical point and it can thus find important appli- 
cations in theoretical astrophysics. Note that fluctuations 
in isothermal spheres close to the critical point have been 
studied by Katz & Okamoto (2000) and Chavanis (2005). 

Finally, we have given an astrophysical application of 
the HMF model in relation with the formation of bars 
in spiral galaxies, following the original idea of Pichon 
and Lynden-Bell. This simple model is consistent with the 
phenomenology of structure formation which results from 
a competition between long-range interactions (gravity) 
and thermal motion (velocity dispersion). It is sometimes 
argued that the moment of inertia of stellar orbits can 
be negative (Pichon 1994). In that case, we must consider 
Eq. with fc < 0. This corresponds to the repulsive 
(anti-ferromagnetic) HMF model. Now, Barre et al. (2002) 
have observed that this model leads to the formation of 
a bicluster. In the stellar disk analogy, the equivalent of 
the "bicluster" would be two bars perpendicular to each 
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other. We do not know whether this type of structure is 
observed in astrophysics. 
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of the eigenvalue close to the 



A Estimate 
critical point 



For simplicity, we restrict ourselves to symmetrical per- 
turbations with respect to the x-axis. The general case is 
treated in Appendix IfI bv another method. In the canon- 
ical ensemble, we have to study the eigenvalue equation 
{V = 0): 



d f 1 



dq 



de\p{e)dej 2nT 



2ir 



q{0') cos{0 ~ e')de' 



2Xq, 
(268) 



where p{9) is the equilibrium density profile and q{9) is 
the perturbation. We consider the clustered phase close 
to the critical point Tc. Thus, we can perform a system- 
atic expansion of the parameters in powers of _B <C 1 or, 
equivalently, in powers of a; = PB <C 1. 

Using Eqs. H20|) and H21f) . the density profile can be 
expanded as 



1 

P 



2tt 
M 



X cos ( 



-(H-2cos^0) 



(269) 



Substituting this result in Eq. (|268|l . and using the expan- 
sion H30|l of the temperature, we obtain 



1 + X cos 6*+— (l + 2cos2 6) 



dq 
d9 



2 

1 + ^ 



q{9') cos{e - e')d9' ^ px^q. (270) 



where we have set A A/ — jix^ with /i = 0(1)- Here, A refers 
to the largest eigenvalue of Eq. H268|l which is equal to 
zero for a; = (the other eigenvalues are A„M = — tto^ -|- 
0{x) for n > 2). Furthermore, the following expansion 
shows that the term of order x vanishes so we have directly 
written A ^ a;^. According to Eq. %\{)\ . we have 



(271) 



where p, has to be determined self-consistently. We thus 
expand the perturbation as 



q{d) = qo{e) + X qi{e) + x" q2{9) + 
and we introduce the differential operator 
d^q 



Cq = TT 



dO^ 



q{e')cos{9~e')de'. 



(272) 



(273) 



To order 0, we have 

Cqo = 0, (274) 
yielding go = sin 6*. To order 1, we get 

= 7rsin(26l), (275) 

yielding 

gi = -isin(26l) + Csin6', (276) 

where C is an arbitrary constant. Finally, to order 2, we 
have after simplification 

Cq2= hi + ^ j sin6l + 7rCsin(26l) - ^sin(36l), (277) 



yielding 



and 



(278) 



C 1 
q2 = Dsine- —sm{2e) + — sin(36'), (279) 

where D is an arbitrary constant. Therefore, close to the 
critical point, the largest eigenvalue of Eq. (|268|l is 

AM = -^x^ or AM = -27r^-|- - 1^. (280) 

We can obtain the expression of the eigenvalue by a 
slightly different method. We consider the Hilbert space of 
27r-periodic continuous real functions with scalar product 



1 



2tt 



(281) 



We note that the operator H273|) is Hermitian in the sense 
that 

{Cf,g)^{f,C9)- (282) 
The equation obtained to second order can be written 

Cq2^gie), (283) 

with 

g{0)= sin6'-|-7rCsin(26') - ^sin(36'). (284) 

We note qo the Kernel of C, i.e. Cqo = 0. Then, we have 
the condition of solvability 



(go,.?) = 0. 



(285) 



Indeed, 



(go, .g) - (go, Cq2) = {Cqo, ga) = 0. (286) 



In our case, go = sin 6, so that the condition (go,<?) = 
with yields Eq. (|775|l . 
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In the microcanonical ensemble, we have to study the 
eigenvalue equation {V ^ 0): 



1 dq 



de \p{e) de J 2t:T 



k 



2-iT 



q{e')cos{e-0')d0' 



2V d$ 

MT'^de 



2\'q. (287) 



To leading order in B ^ 1, the term ^gyr^ can be writ- 
ten ^x'^s'md. Therefore, Eq. I287|l becomes 



d9 



1 + xcose* H (1 + 2cos^ ( 

4 



dq 



+ ^1 + y ^ ^ q{0') cos{e - e')de' = irx^ sin 9 + n'x^q. 

(288) 

The eigenvalue is now 



yielding 



AM = 

4 



^ =-4-" = -T' 



or XM = -IOtt 



- 1 



(289) 



(290) 



Although the onset of instability is the same in the two 
ensembles, the eigenvalues differ in the condensed phase. 



B Some useful identities 

For any system described by a distribution function f{6, v) 
we define the density and the pressure by 



(291) 



The kinetic temperature T — p/ p is equal to the velocity 
dispersion of the particles. If the distribution function is 

of the form / = /(e) with e = ^ + (I>{9), we have 



dp 
d9 



d-P 2 d'P fdf 

f ie)^v dv 



de / dv 



vdv 



which returns the condition of hydrostatic balance (|70|) . 
This relation is equivalent to 



P'ip) 



1 



Now, we note that 



dp 
d^ 



dv. 



(293) 



(294) 



This yields the important identity 
p'ip) -1 



P \T-dv 



(295) 



This identity is valid for both homogeneous and inhomo- 
geneous systems. It may be useful to rederive it in the case 
of homogeneous systems since <P = in that case so that 
the above precedure is ill-determined. 

We consider a homogeneous system described by the 
distribution function / = F{(m^ + a) where a is a La- 
grange multiplier taking into account the conservation of 
mass (normalization). In that case, p = p{a) andp = pio:). 
Now, 



dp 
da 



F'{P—+a)v'^ dv 

|;F(4+«)^d. = -^, (296) 



and 



Eliminating a from the foregoing relations, we obtain 
p'ip) -1 



fl^dv' 

J V 



(298) 



which is consistent with Eq. (|295|) . Now, introducing the 
velocity of sound — p'ip) ^-^d using p — M /27r, we get 



M 



1 



2tt f+^rMdv' 

J — OO V 



(299) 



C General dispersion relation for a gaseous 
system 



We consider the Euler equations 
dp 



dt 



V • (pu) = 0, 



P 



- + (u.V)u 



(300) 



(301) 



in D dimensions and for an arbitrary binary potential of 
interaction u(|r — r'|) such that 



^{r,t)^ J u(r-r')p(r',t)d^r'. 



(302) 



We also consider an arbitrary barotropic equation of state 
p — p{p)- Clearly, p = Cst., u = and ^ = is a station- 
ary solution of Eq. I|300(l - H3()2|l provided that / w(x)(i^x = 
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(for the gravitational potential, this is not the case and 
we must advocate the Jeans swindle). We shall restrict 
ourselves to such homogeneous solutions. The linearized 
Euler equations are 



(303) 



(304) 



S^{r,t) ^ J u{r~r')5p{r',t)d°r'. (305) 
These equations can be combined to give 

-clAdp-pAd^^O, (306) 



d'Sp 2 



where we have introduced the velocity of sound = p'{p)- 
We introduce the Fourier transform of the interaction po- 
tential such that 



u(x) = J e'^ ''u{k)d^k. 



(307) 



Taking the Fourier transform of Eqs. (|305|l and (|306|l with 
the convention Sp ~ g*(k r-ijt)^ combining the result- 
ing expressions, we find that 



uj'^ = clk^ + (27r)^?i(k)fc^p 



(308) 



which is the required dispersion relation. The system will 
be unstable if 



ci + (27r)^u(k)p < 0. 



(309) 



In particular, it is necessary that u(k) < correspond- 
ing to attractive potentials. In that case, the condition of 
instability reads 



< {cl)c„t = (27r)^p|u(k)|, 



(310) 



Then, the unstable lengthscales are determined by Eq. 
H3U9|) . Various situations can happen (see Chavanis 2GG4b) 
depending on the form of the potential {t(k). For the grav- 
itational interaction in D — 3, using Au — 47rG(5(x), we 
have (27r)'^{t(k) — —AirG/k'^. We recover the usual formula 
(Binney & Tremaine 1987) 



^2 ^ clk^ _ 47rGp. 



(311) 



The system is always unstable {{c'^)crit — oo) for suffi- 
ciently large wavelengths or equivalently for 



1/2 



(312) 



where kj is the Jeans wave number. For the HMF model, 
■u„ = if n ^ ±1 and u±i = — In that case, Eq. (|308|l 
returns Eq. (|87ll . The system is unstable for the mode 
n = 1 if c2 < fcM/47r. 



D General dispersion relation for a stellar 
system 

We consider the Vlasov equation 



df df df 

— + V- — -V^- — = 0, 

dt dr dv 



(313) 



in D dimensions and for an arbitrary binary potential of 
interaction u(|r — r'|) as before. Clearly, a distribution 
function / = /o(v) which depends only on the velocity 
is a stationary solution of Eq. (|313|l under the same as- 
sumptions as before. We shall restrict ourselves to such 
homogeneous solutions. The linearized Vlasov equation is 

^+v.f^-V^^.^^O, (314) 
ot or ov 

which must be supplemented by 

S^{r,t)= J u{r-r')6p{r',t)d"r'. (315) 

Taking the Fourier transform of Eqs. (|314(l and (|315|l with 
the convention 6f ~ g»(k r-i^t) ^^^^^ combining the resulting 
expressions, we find that 

(Lu-k- w)6f + (27r)^?i(k)k • ^ / 5f(k, v, t^)d^v = 0. 

dv J 

(316) 

This can be rewritten 



e(k,uj) EE 1 -I- (27r)'°u(k) 



k- ^ 

^d^v^O, (317) 

w — k • V 



which is the dispersion relation of the system. For the 
gravitational interaction in D = 3, using (27r)^M(k) = 
— 47rG/fc2, we recover the usual formula (Binney & Tremaine 
1987) 



47rG r k-^ 



dy_jD 



-k- 



d^v = 0. 



(318) 



For the HMF model, we recover Eq. H117|) . 



E General dispersion relation for a Brownian 
system 

We finally consider a gas of Brownian particles in interac- 
tion described in the strong friction limit by the Smolu- 
chowski equation 



dp 
dt 



{TVp + p\7<P) 



(319) 



where 'P is given by Eq. H3U2|) . Under the same conditions 
as in Appendix 121 the linearized Smoluchowski equation 
is 



dSp 
'dt 



= V • 



^{TV5p + pV5$) 



(320) 
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Taking the Fourier transform of Eq. 1)320(1 , we obtain the 
dispersion relation 



(321) 



which can be compared to Eq. (|308|l . In the gravitational 
case, we get 



i^oj = Tk^ - AnGp. 



(322) 



The stability condition coincides with the Jeans criterion 
for an isothermal gas 



k < k] 



1/2 



(323) 



However, the stable modes are exponentially damped in 
the case of Brownian particles while they have an oscilla- 
tory nature in the case of gaseous systems. For a sinusoidal 
potential of interaction, we recover Eq. (|217|l . 



The equation for p is readily solved and we obtain 
2A{T, - T)e^iT.-T)m 



p{t) = 



1 -L 2{T^-T)t/i 

-L I ^rp 



(331) 



(332) 



Substituting this result in Eq. (|33l)|) and solving the re- 
sulting equation, we get 



X{t) = 



B 



(333) 



The modes ai and a_i are deduced from X and p by 
solving the second order equation a? — Xa+p ~ yielding 



a±i{t) 



X{t) ± ^A{f) 



(334) 



F Generalization of the analytical solution 
0242|t to arbitrary perturbations 

In Sec. 16.41 we have restricted our analysis to density pro- 
files that are always even, i.e. p{—9,t) = p{9,t). It is not 
difficult to relax this hypothesis. Let us write the density 
profile in the general form 



+ 00 

P= ^ a„(i)e*' 

71 — — OO 



(324) 



Substituting this decomposition in Eq. (|234ll . we find that 



M 

"°^2^' 



(325) 



+ Tn a„ = -7i(aia„_i - a_ia„+i). (326) 

To first order in T), — T <C T^, we can neglect a„ with 
\n\ > 3. We thus get 



a±2 = ^«±i- 
The equations for the modes a±i are therefore 



(327) 



^^ + iT-T^)a^^=-^a\a,. (329) 

At that point, it is convenient to introduce the variables 
p = aiQ-i and X = ai +a-i. They satisfy the differential 
equations 



8T 



pX, 



(330) 



The discriminant can be written 



A = 



-8\A\{T,-T) 
|^ + e-2(Tc-T)t/?|- 



(335) 



The constants of integration A and B, which can be pos- 
itive or negative, are fixed by the initial condition. They 
must satisfy A{Q) > 0, i.e. B^ > 8\A\{Tc - T). Then, 
A{t) > at all times. If initially a_i(0) — ai(0), then 
Z\(0) = 0. By Eq. (|^ . this iniphes that A{t) = for 
all times. Therefore, if the initial perturbation is even, it 
remains even during all the evolution. We thus recover 
the results of section ing with ai equals to 2a±i with the 
present notations. Finally, for T = T^, Eq. (|Mn|l and p!?T|l 
lead to 



1 

p(0) 



X{t) = 



X{0) 



1 



k^p{0)t 



(336) 



We can also use this approach to study the dynami- 
cal stability of stationary solutions of Eq. H234|) . By trun- 
cation to order 2 (03 — a_3 ~ 0), we obtain the equi- 
librium relations : B = — fcaf (we suppose = a'Li 
without lack of generality), — alj ~ f (^c ^ T) and 
B^ — 8Tc{Tc — T). These expressions are valid up to order 
1 in Tc — T. Let us compute the four smallest eigenval- 
ues corresponding to the relaxation towards equilibrium. 
Defining a„ — afj -I- exp(At)(5a„, we obtain after lineariza- 
tion MY — where F is a row vector of components 
^(Jfli, (5a_i, (5a2, (5a_2) and 



M = 



/ -i\ + AT 
-AT 



-AT 



-^/W^ \ 

iX + AT -^2TcAT 

^^2TcAT - 4T 

-4V2TcZ\r -CA - 4T y 



where we have set AT — T^. — T. The eigenvalues are the 
zeros of the determinant of M . The two lowest eigenvalues 
scale as (Tc — T) so we set A — |-(Tc — T)Ai. First dividing 
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the two first rows and the two first columns by \J (T^ — T) , 
and dividing the two last rows by 2\/2Tc and the two last 
columns by \/2T^ we obtain at leading order 



det 



'-Ai + 1 -1 -1 
-1 -Ai + 1 -1 
2 0-10 
-2 -1, 



which gives + 2Ai = 0. The smallest eigenvalue Ai = 
is the neutral mode associated to the angle rotation invari- 
ance, whereas the smallest non zero eigenvalue is Ai = — 2 
or equivalently A = —^{T,, — T). The two consecutive 
eigenvalues are of order 0. Considering the determinant 
of M at order 0, we obtain two degenerate eigenvalues 
A = — 4Tc/^. We note that this degeneracy will be re- 
moved at next order. 
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